Sexual selection and speciation: a macro-evolutionary approach

Supplementary Material

library(ggplot2)
library(pander)
library(knitr)
library(tidyr)
library(mgcv)
library(diversitree)
library(repmis)
library(lme4)
library(ape)
library(geiger)
library(grid)
library(gridExtra)
library(nlme)
library(phytools)
library(brms)
library(ggridges)
library(caper)
library(purrr)
library(reshape2)
library(ggExtra)
library(grid)
library(dplyr)
library(BAMMtools)
library(gridExtra)
library(car)
library(coda)
library(MuMIn)
library(parallel)
library(HDInterval)
library(phangorn)
library(kableExtra)
library(RColorBrewer)
library(tibble)
library(stargazer)
library(ggraph)
#devtools::install_github("Ax3man/phylopath")
library(phylopath)
library("EBImage") # for images
library("ggtree")
source("bamm_extinction/functions/check_and_fix_ultrametric.R")
source("functions/essim.R")
source("functions/diversification_rate_calculation.R")
source("functions/Vdodge_function.R")
source('functions/symmetrise_scale.R')
source('functions/Correlation_matrix.R')
source('functions/gheatmap2.R') #Gets rid of spaces between tiles

Here we provide the supplementary material for our study on sexual selection and extinction risk in passerine birds.

Compiling the data

For this project we aimed to compile data on a large set of species in the Passerine bird order that would allow us to account for many evolutionary and environmental predictors of speciation. The following details the source and use of the data.

Speciation Rate

To obtain speciation rate estimates we used Bayesian Analysis of Macroevolutionary Mixtures (BAMM). Additionally, we compared our BAMM output files (event data) against a pre-existing extinction rate dataset (Harvey et al. 2017).

Rather than obtaining one estimate of speciation rate from one tree, we ran BAMM 100 times on 100 trees plus an MCC tree to obtain uncertainty estimates of the speciation rate generated from the phylogeny created by Jetz et al. (2012). We also estimated extinction rate from BAMM.

Sexual Selection

Our proxy for sexual selection was based on measures of sexual dichromatism. Notably, we use two existing datasets; a complete set of male and female plumage scores in Passerines from Dale et al. (2015) and reflectance data from 1000 birds from Armenta et al. (2008).

Briefly, Dale et al. (2015) obtained male and female plumage score, which can be divided to obtain a sexual dichromatism index (SDi). The plumage score was scored using the Handbook of the Birds of the World Del Hoyo et al. (2011). By scanning images of males and females in this book across multiple patches Dale et al. (2015) were able to obtain mean plumage scores from RGB values. Notably concern about the accuracy of these measurements were raised, thus this dataset was previously compared with a more thorough reflectance dataset of Australian birds, the two datasets correlated but was expectedly noisy. Here we compare the Dale et al. (2015) data against another reflectance dataset from Armenta et al. (2008) using colour discriminability as the parameter.

Alongside sexual dichromatism data there is also partial data for other metrics analysed in Dale et al. (2015) such as: Body size, tropical life history, Sexual selection, Cooperative breeding and Migration. We use a phylogenetic PCA (PC1) that estimates male-bias sexual selection from a species level of sexual size dimorphism, social polygyny and paternal care.

Environmental predictors

We obtained species range distribution maps from Birdlife International (through IUCN) (BirdLife International and Handbook of the Birds of the World 2017). These species range maps cover nearly all species of birds. From these spatial information files we were able to obtain estimates of the following (note that not all the data extracted was subsequently used in the analysis):

• Range Size
• Average and standard deviation in 19 bioclimatic variables (each range randomly sampled 1000 times)
• Average and standard deviation in 19 bioclimatic variables from the last-glacial maximum (LGM) and the last-inter-glacial (LIG)
• Net primary productivity (NPP) estimates and variability
• Average and standard deviation of human population density in the species ranges.

The code for extracting the environmental data from the species is provided in another R-markdown document Github link, however no raw data is provided alongside this file due to the large file size of the shapefile and raster information. We have put the extracted environmental variables in a csv filed called complete.dataframe.csv for convenience. Briefly for each of the ~ 6,000 species, we:

1) Randomly sampled the range polygon 1000 times:

#Increse the iterations (defaukt is 4) so we can obtain complete samples of each range
bird.points <- lapply(bird.ranges, 
                       function(x) {spsample(x, n=1000, type="random", iter = 30)})

2) Extracted the bioclim or other (e.g. NPP) data from each of the points:

bird.values <- lapply(bird.points, function(x) {raster::extract(bioclim_data, x)})
bird.values <- as.list(data.frame((bird.values)))

3) Summarised the extracted data across the 1,000 points into a summary value of interest (means and se) for each species and exported that data.

#Obtain means
bird.frame <- lapply(bird.values, function(x) {as.data.frame(x)})
bird.summary <- lapply(bird.frame, function(x) {
(as.data.frame(apply(x,2,mean, na.rm =T )))})

#transpose
bird.means <- t(as.data.frame(bird.summary))
bird.means <- (split(bird.means, 1:19))

#Now add column of species name: Same order carried through
bird.means <- cbind.data.frame(shps.jetz, bird.means)

rownames(bird.means) <- NULL
colnames(bird.means) <- c("binomial","bioclim1", "bioclim2", "bioclim3", "bioclim4", "bioclim5", "bioclim6", "bioclim7", "bioclim8", "bioclim9", "bioclim10", "bioclim11", "bioclim12", "bioclim13", "bioclim14", "bioclim15", "bioclim16", "bioclim17", "bioclim18", "bioclim19")

#Write csv
write.csv(bird.means, 'data/bird.means.csv')

4) This process was repeated when we used gridded data for the LGM, LIG, NPP and human population density, while range size was obtained from the following code:

bird.range.size <- sapply(bird.ranges,
                       function(x) {(area(x))})
bird.range.size <- sapply(bird.range.size, function(x) {sum(x)})
bird.range.size <- as.data.frame(bird.range.size)
bird.range.size <- cbind.data.frame(shps.jetz, bird.range.size)
rownames(bird.range.size) <- NULL
colnames(bird.range.size) <- c("binomial", "range.size.m2")

#write.csv
write.csv(bird.range.size, 'data/bird.range.size.csv')

Generating biologically relevent predictors

Here we assess the relationship between sexual selection and extinction risk. However, in doing so we attempt to take into account as many other predictors of extinction as possible, primarily through environmental variables. Basically, our model structure seeks to contain:

Extinction/Diversification ~ Sexual selection
+ Range size
+ Short temporal variability of temperature (mean BIOCLIM4)
+ Spatial variability of temperature (PCA1) [residual.PCA1]
+ Long-term variability of temperature (LIG)
+ NPP

These predictors can be obtained from the compiled datasets seen here:

plumage.scores <- read.csv('data/plumage_scores.csv')
#Generate sexual dichromatism score: 
plumage.scores$SDi <- abs(plumage.scores$Male_plumage_score - plumage.scores$Female_plumage_score)

#Make DF with enviro variables and plumage scores
complete.dataframe <- read.csv('data/complete.dataframe.csv')
complete.dataframe <- left_join(plumage.scores %>% dplyr::select(binomial, Male_plumage_score, Female_plumage_score, SDi), complete.dataframe, by = "binomial")

Spatial Variability PCAs

To obtain estimates of spatial variability we can obtain PCAs of the standard errors of the bioclimatic variables (excluding seasonality in temperature and precipitation). From there we can run GAM models with the PCA components and range size and extract residuals. We do this as we expect variability will increase alongside increased range size, thus taking the residuals allows us to correct for this association.

Table S1: Loadings for the first two PCs, when we conduct a PCA on the variation (standard error, where n ~ 1,000) pf each bioclim variable except the standard error of temperature seasonality (se.bioclim4) and the standard error of precipitation seasonality (se.bioclim15), as these are already measures of variability and based on other bioclim variables. We use the standard error as, while we attempted to sample each species range 1,000 times, some points returned NA values and would thus affect the standard deviation, this inconvenience should not affect standard error.

restricted.data <- complete.dataframe %>% drop_na(bioclim1) #Drop species without environmental data
PCA.bioclim <- prcomp(restricted.data[c('se.bioclim1', 'se.bioclim2', 'se.bioclim3', 'se.bioclim5', 'se.bioclim6', 'se.bioclim7', 'se.bioclim8', 'se.bioclim9', 'se.bioclim10', 'se.bioclim11', 'se.bioclim12', 'se.bioclim13', 'se.bioclim14', 'se.bioclim16', 'se.bioclim17', 'se.bioclim18', 'se.bioclim19')], #Select se.bioclim
                      scale = TRUE, center = TRUE)
PCA.predictions <- predict(PCA.bioclim)
restricted.data <- cbind(restricted.data, PCA.predictions)

#Create table with PC1 and PC2
as.data.frame(PCA.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1", "PC2")) %>% pander()

	PC1	PC2
se.bioclim1	-0.3328	0.08264
se.bioclim2	-0.2383	0.05458
se.bioclim3	-0.2317	-0.06584
se.bioclim5	-0.312	0.08091
se.bioclim6	-0.3288	0.08275
se.bioclim7	-0.2934	0.09748
se.bioclim8	-0.2916	0.1181
se.bioclim9	-0.3187	0.1125
se.bioclim10	-0.3111	0.06704
se.bioclim11	-0.3296	0.08775
se.bioclim12	-0.1386	-0.4246
se.bioclim13	-0.1596	-0.3083
se.bioclim14	-0.01133	-0.3857
se.bioclim16	-0.1708	-0.3254
se.bioclim17	-0.01923	-0.3932
se.bioclim18	-0.1433	-0.3066
se.bioclim19	-0.01797	-0.3815

#run GAM
PC1.gam <- mgcv::gam(PC1*(-1) ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian") #Inverse as the PC1 loads to the negative (COUNTER-NTUITIVE)

#take residuals
restricted.data$residuals.PC1 <- residuals.gam(PC1.gam) 

#We can do the same for PC2
PC2.gam <- mgcv::gam(PC2 ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian")

#Take residuals
restricted.data$residuals.PC2 <- residuals.gam(PC2.gam)

#plot
par(mfrow = c(1,2))
plot.gam(PC1.gam, residuals = T, main = 'PC1 residuls (temp)')
plot.gam(PC2.gam, residuals = T, main = 'PC2 residuals (precip)')

Figure S1: The relationship between spatial variability in the temperature components (PC1) and log-range size is relatively strong. But not as strong for spatial variability in precipitation (PC2). In the analysis we only used the residuals from PC1

Long term climate variability (LIG anomaly)

To gain estimates for change in climate over the past ~130,000 years we can use the difference in bioclim variables between the LIG and present values. The plots show the two PCs, broadly representing temperature and prescipitation. Data has also been provided for the last-glacial maximum (LGM) but was not used in the analysis.

Table S2: Loadings for the first two PCs of each of the PCAs for the difference in bioclimatic variables between today and the LIG. Here, PC1 is more heavily loaded for absoluted temperature difference, while PC2 is more heavily loaded for absolute difference in precipitation.

historical.variation.data <- as.data.frame(restricted.data[1])

#FOR LIG

historical.variation.data$bio1.LIG.diff <- abs(restricted.data$bioclim1 - restricted.data$LIG.bi1)
historical.variation.data$bio2.LIG.diff <- abs(restricted.data$bioclim2 - restricted.data$LIG.bi2)
historical.variation.data$bio3.LIG.diff <- abs(restricted.data$bioclim3 - restricted.data$LIG.bi3)
historical.variation.data$bio4.LIG.diff <- abs(restricted.data$bioclim4 - restricted.data$LIG.bi4)
historical.variation.data$bio5.LIG.diff <- abs(restricted.data$bioclim5 - restricted.data$LIG.bi5)
historical.variation.data$bio6.LIG.diff <- abs(restricted.data$bioclim6 - restricted.data$LIG.bi6)
historical.variation.data$bio7.LIG.diff <- abs(restricted.data$bioclim7 - restricted.data$LIG.bi7)
historical.variation.data$bio8.LIG.diff <- abs(restricted.data$bioclim8 - restricted.data$LIG.bi8)
historical.variation.data$bio9.LIG.diff <- abs(restricted.data$bioclim9 - restricted.data$LIG.bi9)
historical.variation.data$bio10.LIG.diff <- abs(restricted.data$bioclim10 - restricted.data$LIG.bi10)
historical.variation.data$bio11.LIG.diff <- abs(restricted.data$bioclim11 - restricted.data$LIG.bi11)
historical.variation.data$bio12.LIG.diff <- abs(restricted.data$bioclim12 - restricted.data$LIG.bi12)
historical.variation.data$bio13.LIG.diff <- abs(restricted.data$bioclim13 - restricted.data$LIG.bi13)
historical.variation.data$bio14.LIG.diff <- abs(restricted.data$bioclim14 - restricted.data$LIG.bi14)
historical.variation.data$bio15.LIG.diff <- abs(restricted.data$bioclim15 - restricted.data$LIG.bi15)
historical.variation.data$bio16.LIG.diff <- abs(restricted.data$bioclim16 - restricted.data$LIG.bi16)
historical.variation.data$bio17.LIG.diff <- abs(restricted.data$bioclim17 - restricted.data$LIG.bi17)
historical.variation.data$bio18.LIG.diff <- abs(restricted.data$bioclim18 - restricted.data$LIG.bi18)
historical.variation.data$bio19.LIG.diff <- abs(restricted.data$bioclim19 - restricted.data$LIG.bi19)

historical.variation.data <- historical.variation.data %>% drop_na(bio1.LIG.diff)

#Run a PCA of the difference removing the variation variables (4 and 15)

#For LIG
PCA.LIG.bioclim <- prcomp(historical.variation.data[c(
'bio1.LIG.diff',
'bio2.LIG.diff',
'bio3.LIG.diff',
'bio5.LIG.diff',
'bio6.LIG.diff',
'bio7.LIG.diff',
'bio8.LIG.diff',
'bio9.LIG.diff',
'bio10.LIG.diff',
'bio11.LIG.diff',
'bio12.LIG.diff',
'bio13.LIG.diff',
'bio14.LIG.diff',
'bio16.LIG.diff',
'bio17.LIG.diff',
'bio18.LIG.diff',
'bio19.LIG.diff')],
                      scale = TRUE, center = TRUE)

#Create table with PC1 and PC2
as.data.frame(PCA.LIG.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1.LIG", "PC2.LIG")) %>% pander()

	PC1.LIG	PC2.LIG
bio1.LIG.diff	-0.2956	-0.01943
bio2.LIG.diff	-0.2442	0.0745
bio3.LIG.diff	-0.1436	0.2085
bio5.LIG.diff	-0.2867	-0.2566
bio6.LIG.diff	-0.4094	0.03346
bio7.LIG.diff	-0.4012	-0.09706
bio8.LIG.diff	0.07108	-0.2385
bio9.LIG.diff	-0.3365	-0.00753
bio10.LIG.diff	-0.07014	-0.3766
bio11.LIG.diff	-0.4201	-0.004
bio12.LIG.diff	-0.05181	0.3337
bio13.LIG.diff	-0.2552	0.2591
bio14.LIG.diff	0.05701	0.3298
bio16.LIG.diff	-0.1946	0.3164
bio17.LIG.diff	0.03896	0.3603
bio18.LIG.diff	0.06173	0.277
bio19.LIG.diff	0.08743	0.2856

#Now we can predict the PCA results: 
PCA.LIG.predictions <- as.data.frame(predict(PCA.LIG.bioclim)*(-1)) #So that higher numbers mean more variation
PCA.LIG.predictions <- rename(PCA.LIG.predictions, PC1.LIG = PC1,
       PC2.LIG = PC2)

#Bind them to dataframe, taking only the first two PCAs
historical.variation.data <- cbind(historical.variation.data, PCA.LIG.predictions[1:2])

#Now to the restricted dataframe
restricted.data <- right_join(restricted.data, historical.variation.data %>% dplyr::select(binomial, PC1.LIG, PC2.LIG), by = 'binomial')

PCA.plot <- historical.variation.data %>% ggplot(aes(x = PC1.LIG, y= PC2.LIG))+
  geom_point(shape = 21)+
  theme_minimal()

PCA.plot.m <- ggExtra::ggMarginal(
  p = PCA.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

grid.arrange(PCA.plot.m)

Figure S2: The distribution of the two PCA compontents for the absolute difference in bioclimatic variables between today and the LIG.

Correlations between environmental predictors

We can check the correlations between environmental predictors that we will use in our model. Specifically we want to know whether the following environmental variables are correlated:

• Range size
  
• Short temporal variability of temperature (mean BIOCLIM4)
  
• Spatial variability of temperature (PCA1) [residual.PCA1]
  
• Long-term variability of temperature (LIG)
  
• NPP

To check some of the correlations between the environmental predictors to be used in the PGLS models we can inspect the following correlation plot with the correlation value plotted:

#Draw plot
restricted.data$log.range.size <- log(restricted.data$range.size.m2)
pairs(restricted.data[,c("log.range.size", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP")], lower.panel=panel.smooth, upper.panel=panel.cor)

Figure S3: The highest correlation is between Long term temperature variation and seasonal variation. However this correlation is not excessively high (0.53) so we can be satisfied in including these predictors within our model.

Sexual Dimorphism

The sexual dimorphism dataset looks to have an overdispersed distribution (see phylogenic plot in manuscript). Unfortunately, transformations do not greatly improve the distribution. Our model fit is expected to be reduced by this but further transformations are not obvious and would risk problems of interperetation.

---

Analysis

Simple speciation measures: Diversification rate (DR) and Node Density (ND)

Diversification rates can be simply estimated from thew ‘tippiness of the phylogeny’. Here a log equal splits (also called DR) value can be achieved that is weighted to more recent diversification. It was used by Jetz et al. (2012) and recently (in comparison with BAMM) in a study ivestigation bird diversity alongside elevation (Quintero and Jetz 2018). We can look at speciation rate using ES-Sim adapted by Harvey Michael et al. (2017) and found here https://github.com/mgharvey/ES-sim.

The functions that generate these tip rates estimates are: calculate.log.es, calculate.nd and calculate.tb. Here they are run on 100 trees plus a full MCC tree of passerine birds based on 2500 samples of the posterior (Jetz et al. 2012). This was done using the phanghorn (Schliep 2010).

#Obtain DR (log(es)) estimates, by calling the calculate.dr function
passerine.trees <- read.nexus('bamm_extinction/data/trees/Passerine_Trees_Full.nex')

#DR, ND and TB can be calculated on the 100 trees by: 
  
# es.list <- sapply(passerine.trees, calculate.log.es)
# nd.list <- sapply(passerine.trees, calculate.nd)
# tb.list <- sapply(passerine.trees, calculate.tb)

#To avoid re-running the time-consuming rates, we can save them and load them: 
es.list <- readRDS("data/es.list.rds")
nd.list <- readRDS("data/nd.list.rds")
tb.list <- readRDS("data/tb.list.rds")

#Also load in our MCC tree
MCC.passerine <- read.tree('data/MCC.passerine.tre') #Tree based on 2500 samples of posterior

To gain an estimate of variation between trees we can take the values and obtain a mean, variance and coeeficient of variation (CV). We obtain a CV value to account for unequal variation at higher estimates. The CV is expressed as the ratio of the standard deviation (\(\sigma\)) to the mean (\(\mu\)) (converting to a percentage not needed):

\(C_V = \frac{\sigma}{\mu} \times 100\)

Furthermore, we can account for overdispersion by obtaining a transformed version of CV. In such a case we can calculate CV for a log-normal distribution:

\({\widehat {c_{\rm {v}}}}_{\rm {ln}}={\sqrt {\mathrm {e} ^{s_{\rm {ln}}^{2}}-1}}\)

Given that our ES value is log-transformed, the calculation of \({\widehat {c_{\rm {v}}}}_{\rm {ln}}\) is perhaps more ideal. As for ND we can see from the a qqnorm plot that value of ND is normally distributed, therefore the regular CV is more appropriate for this distribution.

#For ND
#transpose the list so that the elements are the species
t.nd.list <- nd.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
nd.values <- as.data.frame(t(sapply(t.nd.list, function(x) sapply(x, function(x) (x)))))
nd.values <- nd.values %>% tibble::rownames_to_column("binomial")
nd.summary <- melt(nd.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.nd = mean(value), #use normal mean, not a rate
  var.nd = var(value), #Don't really need this if we use CV
  CV.nd = (sd(value)/mean(value)) #Given normal distribution of ND, standard CV formula is appropriate
)

#For ES
#transpose the list so that the elements are the species
t.es.list <- es.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
es.values <- as.data.frame(t(sapply(t.es.list, function(x) sapply(x, function(x) (x)))))
es.values <- es.values %>% tibble::rownames_to_column("binomial")
es.summary <- melt(es.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.loges = log(1/mean(1/exp(value))), #use harmonic mean for rates and log-transform post calculation of HM
  var.loges = var(value), #Don't really need this if we use CV
  CV.loges = sqrt(exp(var(value)) - 1) #CV for log-normal distribution of ES
)

#For TB
#transpose the list so that the elements are the species
t.tb.list <- tb.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
tb.values <- as.data.frame(t(sapply(t.tb.list, function(x) sapply(x, function(x) (x)))))
tb.values <- tb.values %>% tibble::rownames_to_column("binomial")
tb.summary <- melt(tb.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.tb = mean(value), #use normal mean; not a rate
  var.tb = var(value), #Don't really need this if we use CV
  CV.tb = (sd(value)/mean(value)) #Given normal distribution of tb, standard CV formula is appropriate
)

MCC.nd.df <- as.data.frame(calculate.nd(MCC.passerine)) %>% tibble::rownames_to_column("binomial")
MCC.dr.df <- as.data.frame(calculate.log.es(MCC.passerine)) %>% tibble::rownames_to_column("binomial")

#Join the dataframes
dr.summary <- full_join(nd.summary, es.summary, by = "binomial")
dr.summary <- full_join(dr.summary, tb.summary, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.nd.df, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.dr.df, by = "binomial")

dr.summary <- dr.summary %>% rename(TipLabel = binomial, MCC.ND = nd, MCC.DR = es)

saveRDS(dr.summary, 'data/dr.summary.rds')

dr.summary <- readRDS('data/dr.summary.rds')
#Plot a summary of logES and TB with the DRs with their weightings
dr.plot <- dr.summary %>% ggplot(aes(x = mean.loges, y = mean.nd, 
                          fill = (1/(CV.nd))/sum((1/(CV.nd)))*100, 
                      size = (1/(CV.loges))/sum((1/(CV.loges)))*100))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.loges - 0.674*sqrt(var.loges), xmax = mean.loges + 0.674*sqrt(var.loges)), 
                 size = 0.015)+
  geom_errorbar(aes(ymin = mean.nd - 0.674*sqrt(var.nd), ymax = mean.nd + 0.674*sqrt(var.nd)), 
                size = 0.015)+
  theme(legend.position="bottom")+
  labs(size = 'Weight (%) [using DR]', y='Node Density [ND]', x= 'Log Diversification Rate [DR]', fill = 'Weight (%) [using ND] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "RdPu", direction = 1, breaks=c(0,0.015,.03),
                           limits=c(0,0.03))
marginal.dr <-ggMarginal(
  p = dr.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = '#0000009C',
  fill = '#D12E6769'
)

grid.newpage()
grid.draw(marginal.dr)

Figure S4: \(\lambda_{DR}\) and \(\lambda_{ND}\) are correlated in their estimation of species rate. Furthermore, we see similar variation across each species in the estimates across 100 trees as the inverse of the coefficeint of variation (weight) is approximately equivalent across all values of ND and DR (colours and size). Notably, these estimates speciation rate are relatively uncertain across all species, with 50 % CIs in both axes plotted.

Combine the estimates of DR and ND with the restricted dataframe containing all the other variables. We can inspect the trend in SDi in relation to the speciation rates.

restricted.data <- left_join(restricted.data, dr.summary, by = 'TipLabel')

es.plot <- restricted.data %>% ggplot(aes(y = mean.loges, x = SDi, fill = mean.loges))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "YlOrBr", direction = 1, guide = FALSE)+
ylab("mean.DR")+
theme_minimal()

nd.plot <- restricted.data %>% ggplot(aes(y = mean.nd, x = SDi, fill = mean.nd))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "Greens", direction = 1, guide = FALSE)+
theme_minimal()

tb.plot <- restricted.data %>% ggplot(aes(y = mean.tb, x = SDi, fill = mean.tb))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "PuBu", direction = 1, guide = FALSE)+
theme_minimal()

grid.arrange(es.plot, nd.plot, tb.plot, nrow = 3)

Figure S5: Scatter plots showing the raw relationship between sexual dichromatism (SDi) and speciation rates. Across all three measures the pattern and spread is similar, with no obvious relationship. Note that terminal branch lenght (mean.tb) was not used in the analysis.

BAMM measures of speciation and extinction

Set up BAMM parameters

For the use of BAMM we used the following code to generate the parameters across the 100 trees. Each paramaeter value is specified in this code chunk to ensure reproducibility when re-running on any of the trees. These same parameters were also used for the MCC tree (with the seed set at 2500).

name.passerine.tree <- names(passerine.trees)

priors <- sapply(name.passerine.tree, function(x) {
  setBAMMpriors(passerine.trees[[x]], outfile = NULL)
})

sapply(name.passerine.tree, function(x) {
  write.tree(passerine.trees[[x]], paste("data/bamm_files/", x, ".tre", sep=""))
})

# Here is a block of parameters for the control file. We can make a control file for each tree:
params <- list()
for (x in name.passerine.tree) {

# GENERAL SETUP AND DATA INPUT

params[[x]] <- list(modeltype = 'speciationextinction',
# Specify "speciationextinction" or "trait" analysis
                                  
treefile = paste(x, ".tre", sep=""),
# File name of the phylogenetic tree to be analyzed

runInfoFilename = 'run_info.txt',
# File name to output general information about this run

sampleFromPriorOnly = 0,
# Whether to perform analysis sampling from prior only (no likelihoods computed)

runMCMC = 1,
# Whether to perform the MCMC simulation. If runMCMC = 0, the program will only
# check whether the data file can be read and the initial likelihood computed

loadEventData = 0,                       
# Whether to load a previous event data file

eventDataInfile = 'event_data_in.txt',
# File name of the event data file to load, used only if loadEventData = 1

initializeModel = 1,
# Whether to initialize (but not run) the MCMC. If initializeModel = 0, the
# program will only ensure that the data files (e.g., treefile) can be read

useGlobalSamplingProbability = 1,
# Whether to use a "global" sampling probability. If False (0), expects a file
# name for species-specific sampling probabilities (see sampleProbsFilename)
                                        
globalSamplingFraction = 1,
# The sampling probability. If useGlobalSamplingProbability = 0, this is ignored
# and BAMM looks for a file name with species-specific sampling fractions

sampleProbsFilename = 'sample_probs.txt',
# File name containing species-specific sampling fractions

seed = as.numeric(gsub("tree_", "", x, perl = TRUE)),
# Seed for the random number generator. Set for reproducibility to the number of the treefile

overwrite = 1,
# If True (1), the program will overwrite any output files in the current
# directory (if present)


# PRIORS

expectedNumberOfShifts = 100,
# prior on the number of shifts in diversification
# Suggested values: 
#     expectedNumberOfShifts = 1.0 for small trees (< 500 tips)
#  expectedNumberOfShifts = 10 or even 50 for large trees (> 5000 tips) 
 
lambdaInitPrior = as.numeric(priors['lambdaInitPrior', x]),
# Prior (rate parameter of exponential) on the initial lambda value for rate
# regimes

lambdaShiftPrior = 0.05,
# Prior (std dev of normal) on lambda shift parameter for rate regimes
# You cannot adjust the mean of this distribution (fixed at zero, which is
# equal to a constant rate diversification process)

muInitPrior = as.numeric(priors['muInitPrior', x]),
# Prior (rate parameter of exponential) on extinction rates  

lambdaIsTimeVariablePrior = 1,
# Prior (probability) of the time mode being time-variable (vs. time-constant)
            

# MCMC SIMULATION SETTINGS & OUTPUT OPTIONS

numberOfGenerations = '100000000',
# Number of generations to perform MCMC simulation

mcmcOutfile = 'mcmc_out.txt',
# File name for the MCMC output, which only includes summary information about
# MCMC simulation (e.g., log-likelihoods, log-prior, number of processes)

mcmcWriteFreq = 1000,
# Frequency in which to write the MCMC output to a file

eventDataOutfile = 'event_data.txt',
# The raw event data (these are the main results). ALL of the results are
# contained in this file, and all branch-specific speciation rates, shift
# positions, marginal distributions etc can be reconstructed from this output.
# See R package BAMMtools for working with this output

eventDataWriteFreq = 1000,
# Frequency in which to write the event data to a file

printFreq = 10000,
# Frequency in which to print MCMC status to the screen

acceptanceResetFreq = 1000,
# Frequency in which to reset the acceptance rate calculation
# The acceptance rate is output to both the MCMC data file and the screen

outName = x,
# Optional name that will be prefixed on all output files (separated with "_")
# If commented out, no prefix will be used


# OPERATORS: MCMC SCALING OPERATORS

updateLambdaInitScale = 2,
# Scale parameter for updating the initial speciation rate for each process

updateLambdaShiftScale = 0.1,
# Scale parameter for the exponential change parameter for speciation

updateMuInitScale = 2,
# Scale parameter for updating initial extinction rate for each process

updateEventLocationScale = 0.1,
# Scale parameter for updating LOCAL moves of events on the tree
# This defines the width of the sliding window proposal
 
updateEventRateScale = 4,
# Scale parameter (proportional shrinking/expanding) for updating
# the rate parameter of the Poisson process

# OPERATORS: MCMC MOVE FREQUENCIES

updateRateEventNumber = 1,
# Relative frequency of MCMC moves that change the number of events

updateRateEventPosition = 0.25,
# Relative frequency of MCMC moves that change the location of an event on the
# tree

updateRateEventRate = 1,
# Relative frequency of MCMC moves that change the rate at which events occur 

updateRateLambda0 = 1,
# Relative frequency of MCMC moves that change the initial speciation rate
# associated with an event

updateRateLambdaShift = 1,
# Relative frequency of MCMC moves that change the exponential shift parameter
# of the speciation rate associated with an event

updateRateMu0 = 1,
# Relative frequency of MCMC moves that change the extinction rate for a given
# event

updateRateLambdaTimeMode = 0,
# Relative frequency of MCMC moves that flip the time mode
# (time-constant <=> time-variable)

localGlobalMoveRatio = 10,
# Ratio of local to global moves of events 


# INITIAL PARAMETER VALUES

lambdaInit0 = 0.032,
# Initial speciation rate (at the root of the tree)

lambdaShift0 = 0,
# Initial shift parameter for the root process

muInit0 = 0.005,
# Initial value of extinction (at the root)

initialNumberEvents = 0,
# Initial number of non-root processes


# METROPOLIS COUPLED MCMC

numberOfChains = 1,
# Number of Markov chains to run

deltaT = 0.01,
# Temperature increment parameter. This value should be > 0
# The temperature for the i-th chain is computed as 1 / [1 + deltaT * (i - 1)]

swapPeriod = 1000,
# Number of generations in which to propose a chain swap

chainSwapFileName = 'chain_swap.txt',
# File name in which to output data about each chain swap proposal.
# The format of each line is [generation],[rank_1],[rank_2],[swap_accepted]
# where [generation] is the generation in which the swap proposal was made,
# [rank_1] and [rank_2] are the chains that were chosen, and [swap_accepted] is
# whether the swap was made. The cold chain has a rank of 1.


# NUMERICAL AND OTHER PARAMETERS

minCladeSizeForShift = 3,
# Allows you to constrain location of possible rate-change events to occur
# only on branches with at least this many descendant tips. A value of 1
# allows shifts to occur on all branches. 

segLength = 0.025,
# Controls the "grain" of the likelihood calculations. Approximates the
# continuous-time change in diversification rates by breaking each branch into
# a constant-rate diversification segments, with each segment given a length
# determined by segLength. segLength is in units of the root-to-tip distance of
# the tree. So, if the segLength parameter is 0.01, and the crown age of your
# tree is 50, the "step size" of the constant rate approximation will be 0.5.
# If the value is greater than the branch length (e.g., you have a branch of
# length < 0.5 in the preceding example) BAMM will not break the branch into
# segments but use the mean rate across the entire branch.

outName = x)
  }

bammcontrolfile <- list()
for (x in name.passerine.tree) {
  bammcontrolfile[x] <- paste("data/bamm_files/control_", x, ".txt", sep="")
}

# Now writing control parameters to file

for (x in name.passerine.tree) {generateControlFile(file = bammcontrolfile[[x]], type = "diversification", params = params[[x]])}

Run analysis

BAMM can be run through the terminal through the following syntax: bamm -c control_tree_xxxx.txt. To generate these commands we can use a loop function, from which we get:

bamm.commands <- list()
for (x in name.passerine.tree) {
  bamm.commands[x] <- paste("bamm -c control_", x, ".txt", sep="")
}

The analysis was run over multiple CPU’s, each generating a respective MCMC and EventData output. Due to the size of the event data file (~ 50 Gb in total) they are not included as supplementary material here. However we have simplified the event data objects into tip rate estimates of the mean and variance across 100 trees + MCC.

Read in the event data and extract tip data

Now we have a series of trees and event data we can read in. Firstly and most importantly let’s check for convergence in the MCC tree through the following:

Table S3: Effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) for the run on the MCC indicate that BAMM converges (ESS > 200).

#Read in the tree and MCMC to check for convergence 
MCC.BAMM.tree  <- read.tree("data/BAMM_MCC/MCC.passerine.tre") #Same as other MCC tree already loaded
MCC.BAMM.mcmc <- read.csv( "data/BAMM_MCC/tree_MCC_mcmc_out.txt" , stringsAsFactors=F)

#Plot of convergence can be generated by:
#plot(MCC.BAMM.mcmc$logLik ~ MCC.BAMM.mcmc$generation)

#Looks like it has converged so let's discard burn in: 
burnstart <- floor(0.1 * nrow(MCC.BAMM.mcmc))
postburn <- MCC.BAMM.mcmc[burnstart:nrow(MCC.BAMM.mcmc), ]

#We can also check effective population sizes of the log-likelihhod and number of shift events in each sample
#We want at least 200 (although that's on the low side)

cbind(effectiveSize(postburn$N_shifts), effectiveSize(postburn$logLik)) %>% `colnames<-`(c("N_Shifts", "logLik")) %>% `rownames<-`("Effective Sample Size") %>% pander()

	N_Shifts	logLik
Effective Sample Size	2456	1302

Those values are well above 200, so we are good to use the MCC data in analysis! We can also check the convergence of BAMM across 100 runs of BAMM. The Raw MCMCs are not included in this file or repository but we can read in a dataframe that has extracted the ESS for each of the runs.

Table S4: For the 100 trees that BAMM was run on effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) also indicates that BAMM converges with the minimum for each parameter across the 100 trees being over 200.

ESS <- readRDS('data/ESS.rds')
summary(ESS) %>% pander()

N_Shifts	logLik
Min. : 723.4	Min. : 278.7
1st Qu.:1396.1	1st Qu.: 511.3
Median :1618.6	Median : 638.6
Mean :1709.5	Mean : 669.9
3rd Qu.:2070.8	3rd Qu.: 804.8
Max. :3425.6	Max. :1244.1

Given that it appears BAMM converges we can we can make a data frame with the mean and variance for extinction and speciation tip rates from the large event data set for the MCC with the following code and then plot the variation we see in tip-rates.

# Read in Event Data
MCC.BAMM.ED <- getEventData(MCC.BAMM.tree,  "data/BAMM_MCC/tree_MCC_event_data.txt", burnin=0.1, nsamples=1000)

## Reading event datafile:  data/BAMM_MCC/tree_MCC_event_data.txt 
##      ...........
## Read a total of 100000 samples from posterior
## 
## Discarded as burnin: GENERATIONS <  9999000
## Analyzing  1000  samples from posterior
## 
## Setting recursive sequence on tree...
## 
## Done with recursive sequence

saveRDS(MCC.BAMM.ED, 'data/MCC.BAMM.ED.rds')
#From the Event Data we can extract 

library(purrr)
#BAMM.EventData <- readRDS('data/BAMM.EventData.rds')

#Big lapply over each tree in BAMM event data
BAMM.extraction.function <- function(x) {
######Get mean and var for lambda
#Transpose list so each element in the list is a species
transposed.lambda <- lapply(purrr::transpose(x$tipLambda), unlist)

#Now turn it into a df with mean and variance
lambda <- sapply(transposed.lambda, function(x) {
  mean.lambda = mean(log(x))
  var.lambda = var(log(x))
  return(c(mean.lambda, var.lambda))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.lambda", "var.lambda"))

lambda$TipLabel <- x[["tip.label"]]

#####NOW FOR Extinction

#Transpose list so each element in the list is a species
transposed.mu <- lapply(purrr::transpose(x$tipMu), unlist)

#Now turn it into a df with mean and variance
mu <- sapply(transposed.mu, function(x) {
  mean.mu = mean(log(x))
  var.mu = var(log(x))
  return(c(mean.mu, var.mu))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.mu", "var.mu"))

mu$TipLabel <- x[["tip.label"]]

left_join(lambda, mu, by = "TipLabel")
}

MCC.BAMM.df <- BAMM.extraction.function(MCC.BAMM.ED)

#Save df for later use
saveRDS(MCC.BAMM.df, 'data/MCC.BAMM.df.rds')

MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#Plot a summary of logES and TB with the DRs with their weightings
BAMM.MCC.plot <- MCC.BAMM.df %>% ggplot(aes(x = mean.lambda, y = mean.mu, 
                          fill = 1/var.lambda, 
                      size = 1/var.mu))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.lambda - 0.674*sqrt(var.lambda), xmax = mean.lambda + 0.674*sqrt(var.lambda)),
                 size = 0.0025)+
  geom_errorbar(aes(ymin = mean.mu - 0.674*sqrt(var.mu), ymax = mean.mu + 0.674*sqrt(var.mu)),
                size = 0.0025)+
  # scale_y_continuous(trans = "log")+
  # scale_x_continuous(trans = "log")+
  # xlim(-8, 3)+
  # ylim(-10,2)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Inverse log(var) / Weight [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Inverse log(var) / Weight [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1)

BAMM.variance <- ggExtra::ggMarginal(
  p = BAMM.MCC.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = '#BA3B1C91'
)

grid.newpage()
grid.draw(BAMM.variance)

Figure S5: The tip-rate estimate for BAMM are highly variable within each run of BAMM. Across most species there is high variability in the posterior distribution of tip-rate estimates. Here we show mean values, weights based on the variance and 50 % CIs.

Analysis of BAMM results

Given the high variability in tip-rate estimates from the above plot, below we peform some diagnostics on BAMM to demonstrate that the variability is unlikely an error in sampling or parameters, rather an inherent aspect of BAMM. The following is not inherently necessary to understandfor the conclusions drawn from the paper, however they do raise a set of methodological questions about BAMM that warrant further investigation.

Credible number of shifts

To plot the credible shift set, we need the prior distribution on the number of rate shifts (this is generated internally by BAMMtools). We can then estimate the credible set of rate shifts using the BAMMtools function credibleShiftSet:

css <- credibleShiftSet(MCC.BAMM.ED, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

#Now we obtain the number of distinct shifts: (Out of 1000 samples this is super high, essentially each one is distinct)

css$number.distinct

## [1] 950

This value indicates that from a posterior sample of 1000 we have 950 unique rate shifts, which is very high. Arguably, across such a large phylogenetic tree, even if the number of shifts reaches convergence the location of those > 50 shifts will be quite variable, with the variation likely impacting tip-rate estimates. Even if we increase the sampling of the posterior to 2,500 instead of 1,000 we still see that the CSS is high (close to 2500). We can have a look what our top few shifts are:

summary(css)

## 
##  95 % credible set of rate shift configurations sampled with BAMM
## 
## Distinct shift configurations in credible set:  950
## 
## Frequency of 9 shift configurations with highest posterior probability:
## 
## 
##    rank     probability cumulative  Core_shifts
##          1      0.001      0.001         54
##          2      0.001      0.002         44
##          3      0.001      0.003         46
##          4      0.001      0.004         47
##          5      0.001      0.005         49
##          6      0.001      0.006         45
##          7      0.001      0.007         47
##          8      0.001      0.008         43
##          9      0.001      0.009         45
## 
## ...omitted 941 additional distinct shift configurations
## from the credible set. You can access the full set from your 
## credibleshiftset object

Notably, each of the top shifts has low probability, indicating that out of the entire posterior sample we cannot differentiate which shift configuration is more likely. Based on the follwing eastimates of the BayesFactor, we are convinced that the number of shifts is non-zero, however it is still quite a wide distribution, sourced from potentially weakly non-informative priors. But, this is at odds with the ESS and the following plot:

round(computeBayesFactors(MCC.BAMM.mcmc, expectedNumberOfShifts=100, burnin=0.1)[,1], digits = 2)

##      42      43      44      45      46      47      48      49      50 
##    1.00   17.17   24.48   49.45  120.71  283.77  454.33  784.80 1291.85 
##      51      52      53      54      55      56      57      58      59 
## 1813.33 2561.62 3608.07 4574.91 5492.44 6322.11 7141.12 7720.26 7886.28 
##      60      61      62      63      64      65      66      67      68 
## 7846.73 7661.83 7116.15 6466.36 5749.34 4811.16 4136.80 3338.17 2674.70 
##      69      70      71      72      73      74      75      76      77 
## 2137.61 1593.48 1218.40  868.01  643.91  452.36  311.07  199.17  135.99 
##      78      79      80      81      82      83      84      85      87 
##   85.85   65.03   40.87   14.74   13.40    9.02   10.63    1.53    1.56

plotPrior(MCC.BAMM.mcmc, expectedNumberOfShifts=100)

Figure S6: The apparant convergence of the number of shifts in the posterior is at odds with the variability seen in the CSS and although there is greater certainty in the number of shifts being within the range above, the position of those shifts remain highly variable.

We can compare our run of BAMM against Harvey et al. (2017) who used BAMM on a genetic-only MCC tree with different parameters

load('data/Hackett_split_eventsample.rda')
css3 <- credibleShiftSet(ed, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)
css3$number.distinct

## [1] 1188

Based on Their posterior sample of 1250, 1188 are unique. Like our run on the MCC, we find high uncertainity in the position of shift-configurations.

Let’s run a model based on genetic data from the Harvey MCC

Harvey.BAMM.df <- BAMM.extraction.function(ed)

Harvey.BAMM.df %>% ggplot(aes(x = mean.lambda, y = mean.mu, 
                          fill = 1/(var.lambda), 
                      size = 1/(var.mu)))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.lambda - 0.674*sqrt(var.lambda), xmax = mean.lambda + 0.674*sqrt(var.lambda)), 
                 size = 0.0025)+
  geom_errorbar(aes(ymin = mean.mu - 0.674*sqrt(var.mu), ymax = mean.mu + 0.674*sqrt(var.mu)), 
                size = 0.0025)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Weight [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Weight [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1)

# BAMM.variance <- ggExtra::ggMarginal(
#   p = BAMM.MCC.plot,
#   type = 'density',
#   margins = 'both',
#   size = 5,
#   colour = 'black',
#   fill = '#BA3B1C91'
# )

Figure S7: Based on the event data from the BAMM run by Harvey et al. (2017) we find that in both our case and theirs the tip-rate estimates for many species are extrememly variable. here we present mean estimates with 50 % CIs.

Time-constant

We also ran a time-constant MCC of BAMM

#Read in the tree and MCMC to check for convergence 
MCC.BAMM.tree  <- read.tree("data/BAMM_MCC/MCC.passerine.tre") #Same as other MCC tree already loaded
MCC.BAMM.Constant.mcmc <- read.csv( "data/BAMM_MCC/TimeConstant/tree_MCC_Constant_mcmc_out.txt" , stringsAsFactors=F)

#Plot of convergence can be generated by:
#plot(MCC.BAMM.mcmc$logLik ~ MCC.BAMM.mcmc$generation)

#Looks like it has converged so let's discard burn in: 
burnstart <- floor(0.1 * nrow(MCC.BAMM.Constant.mcmc))
postburn <- MCC.BAMM.Constant.mcmc[burnstart:nrow(MCC.BAMM.Constant.mcmc), ]

#We can also check effective population sizes of the log-likelihhod and number of shift events in each sample
#We want at least 200 (although that's on the low side)

cbind(effectiveSize(postburn$N_shifts), effectiveSize(postburn$logLik)) %>% `colnames<-`(c("N_Shifts", "logLik")) %>% `rownames<-`("Effective Sample Size") %>% pander()


# Read in Event Data
MCC.BAMM.ED.Constant <- getEventData(MCC.BAMM.tree,  "data/BAMM_MCC/TimeConstant/tree_MCC_Constant_event_data.txt", burnin=0.1, nsamples=1000)
saveRDS(MCC.BAMM.ED.Constant, 'data/MCC.BAMM.ED.Constant.rds')
#From the Event Data we can extract 

css2 <- credibleShiftSet(MCC.BAMM.ED.Constant, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

#Now we obtain the number of distinct shifts: (Out of 1000 samples this is super high, essentially each one is distinct)

css2$number.distinct

plotPrior(MCC.BAMM.Constant.mcmc, expectedNumberOfShifts=100)

MCC.BAMM.df.Constant <- BAMM.extraction.function(MCC.BAMM.ED.Constant)

MCC.BAMM.df.Constant %>% ggplot(aes(x = mean.lambda, y = MCC.BAMM.df$mean.lambda))+
  geom_point()
cor(MCC.BAMM.df.Constant$mean.lambda, MCC.BAMM.df$mean.lambda)

PGLS Models

The method behind running the PGLS models is as follows:

1: Estimate the phylogenetic signal (\(\lambda\)) by running a model without interactions and all six predictor variables. The value of \(\lambda\) obtained here will be fixed in all successive models. The value is fixed for subsequent models as independently estimating it in each case becomes computationally intensive. Based on preliminary analysis even slight differences in lambda lead to large reductions in AIC.

2: Create a global model of six predictor variables plus the fiver interactions between sexual dichromatism and the other variables.

3: Dredge the global model but fix the six independent predictors, hence conducting model selection on the interaction terms.

4: Take the top model in the MCC model and run it on the 100 phylogenetic trees.

5: Repeat this step for DR, ND, BAMM-speciation and BAMM-extinction.

PGLS Models on DR and ND

Using corPagel we can estimate the phylogenetic signal for a model with all predictors (interactions do not appear to affect the estimate of \(\lambda\)):

#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(restricted.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(restricted.data) <- restricted.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.corPagel <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.DR.corPagel, 'data/MCC.DR.corPagel.rds')

#Run a corPagel model to estimate lambda for ND
MCC.ND.corPagel <- gls(MCC.ND ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.ND.corPagel, 'data/MCC.ND.corPagel.rds')

Inspect the \(\lambda\) value, which can then be fixed for successive models.

MCC.DR.corPagel <- readRDS('data/MCC.DR.corPagel.rds')
MCC.ND.corPagel <- readRDS('data/MCC.ND.corPagel.rds')
MCC.DR.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()

• DR lambda: 0.985

MCC.ND.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("ND lambda") %>% pander()

• ND lambda: 0.9996

The lambda is high and not too dissimilar than if we assume Brownian Motion \(\lambda = 1\). However, given the large sample size this difference may have an effect on the results so, to play safe we included it as a fixed value for \(\lambda\) in all successive models.

#Run model for DR
MCC.DR <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Run model for ND
MCC.ND <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP
                      + SDi*log(range.size.m2)
                      + SDi*bioclim4
                      + SDi*residuals.PC1
                      + SDi*PC1.LIG
                      + SDi*NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("restricted.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.ND.model <- pdredge(MCC.ND, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.ND.model, "data/dredged.ND.model.rds")

dredged.DR.model <- pdredge(MCC.DR, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.DR.model, "data/dredged.DR.model.rds")

Table S5: The dredged models both show the top model is one with no interactions, with \(\delta AICc > 4\) in both cases. We can be resonably confident that interactions are unlikely to affect the pattern of speciation we see in passerine birds.

dredged.DR.model <- readRDS("data/dredged.DR.model.rds")
dredged.ND.model <- readRDS("data/dredged.ND.model.rds")
kable(dredged.DR.model, "html", caption = "DR Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

DR Dredge Table

	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	SDi	bioclim4:SDi	log(range.size.m2):SDi	NPP:SDi	PC1.LIG:SDi	residuals.PC1:SDi	df	logLik	AICc	delta	weight
0	-2.879811	2.14e-05	-0.0065780	3e-07	0.0015864	-0.0024295	-0.0012787	NA	NA	NA	NA	NA	8	-3409.547	6835.118	0.00000	9.988652e-01
8	-2.879495	2.14e-05	-0.0065788	3e-07	-0.0001028	-0.0023765	-0.0012942	NA	NA	NA	0.0003487	NA	9	-3416.089	6850.208	15.09006	5.281325e-04
16	-2.879345	2.15e-05	-0.0066034	3e-07	0.0015509	-0.0015982	-0.0013082	NA	NA	NA	NA	-0.0001638	9	-3416.488	6851.007	15.88895	3.542138e-04
2	-2.884345	2.14e-05	-0.0064052	3e-07	0.0015991	-0.0024487	-0.0003035	NA	-0.0000369	NA	NA	NA	9	-3416.842	6851.714	16.59633	2.486904e-04
1	-2.878432	1.57e-05	-0.0065688	3e-07	0.0015872	-0.0024180	-0.0015875	1.2e-06	NA	NA	NA	NA	9	-3421.210	6860.452	25.33361	3.150533e-06
24	-2.878703	2.16e-05	-0.0066190	3e-07	-0.0004658	-0.0010549	-0.0013437	NA	NA	NA	0.0004121	-0.0002586	10	-3422.848	6865.734	30.61631	2.245227e-07
4	-2.884812	2.14e-05	-0.0065580	7e-07	0.0015560	-0.0024009	-0.0003587	NA	NA	-1e-07	NA	NA	9	-3424.307	6866.644	31.52618	1.424567e-07
10	-2.890455	2.13e-05	-0.0061602	3e-07	-0.0001797	-0.0024196	0.0010665	NA	-0.0000894	NA	0.0003710	NA	10	-3423.322	6866.681	31.56307	1.398533e-07
18	-2.885405	2.15e-05	-0.0063727	3e-07	0.0015667	-0.0015933	-0.0000022	NA	-0.0000495	NA	NA	-0.0001699	10	-3423.771	6867.580	32.46202	8.922136e-08
9	-2.879477	2.13e-05	-0.0065786	3e-07	-0.0000996	-0.0023764	-0.0012983	0.0e+00	NA	NA	0.0003480	NA	10	-3427.745	6875.529	40.41070	1.676618e-09
17	-2.876970	1.25e-05	-0.0065992	3e-07	0.0015377	-0.0012411	-0.0018098	1.9e-06	NA	NA	NA	-0.0002306	10	-3428.008	6876.054	40.93602	1.289324e-09
3	-2.889434	1.34e-05	-0.0061257	3e-07	0.0016199	-0.0024624	0.0007743	1.7e-06	-0.0000939	NA	NA	NA	10	-3428.372	6876.782	41.66339	8.962234e-10
12	-2.884227	2.14e-05	-0.0065599	7e-07	-0.0000916	-0.0023507	-0.0004248	NA	NA	-1e-07	0.0003405	NA	10	-3430.871	6881.780	46.66167	7.362996e-11
26	-2.893478	2.15e-05	-0.0060556	3e-07	-0.0006022	-0.0009958	0.0018498	NA	-0.0001211	NA	0.0004478	-0.0002816	11	-3430.028	6882.101	46.98312	6.269768e-11
20	-2.884117	2.15e-05	-0.0065817	7e-07	0.0015263	-0.0016648	-0.0004369	NA	NA	-1e-07	NA	-0.0001454	10	-3431.269	6882.576	47.45753	4.945782e-11
6	-2.895008	2.14e-05	-0.0061863	7e-07	0.0015803	-0.0024391	0.0018191	NA	-0.0000790	-1e-07	NA	NA	10	-3431.545	6883.128	48.01036	3.751366e-11
25	-2.877827	1.81e-05	-0.0066164	3e-07	-0.0003444	-0.0009513	-0.0015348	7.0e-07	NA	NA	0.0003863	-0.0002784	11	-3434.448	6890.942	55.82373	7.542887e-13
11	-2.891727	1.89e-05	-0.0060861	3e-07	-0.0000889	-0.0024252	0.0013324	5.0e-07	-0.0001044	NA	0.0003533	NA	11	-3434.894	6891.834	56.71575	4.828789e-13
5	-2.883424	1.60e-05	-0.0065495	7e-07	0.0015572	-0.0023903	-0.0006670	1.1e-06	NA	-1e-07	NA	NA	10	-3435.976	6891.989	56.87126	4.467550e-13
19	-2.894722	8.10e-06	-0.0058780	3e-07	0.0015805	-0.0010549	0.0020214	2.8e-06	-0.0001542	NA	NA	-0.0002814	11	-3435.056	6892.158	57.03960	4.106919e-13
28	-2.882980	2.16e-05	-0.0065992	6e-07	-0.0004288	-0.0011303	-0.0005650	NA	NA	-1e-07	0.0004000	-0.0002393	11	-3437.666	6897.377	62.25931	3.020435e-14
14	-2.901190	2.13e-05	-0.0059396	7e-07	-0.0002030	-0.0024099	0.0032037	NA	-0.0001318	-1e-07	0.0003719	NA	11	-3438.023	6898.091	62.97260	2.114363e-14
22	-2.895506	2.14e-05	-0.0061663	7e-07	0.0015517	-0.0016637	0.0019999	NA	-0.0000885	-1e-07	NA	-0.0001541	11	-3438.494	6899.034	63.91561	1.319496e-14
27	-2.898432	1.32e-05	-0.0057816	3e-07	-0.0003798	-0.0007224	0.0029159	1.7e-06	-0.0001787	NA	0.0004037	-0.0003394	12	-3441.451	6906.955	71.83674	2.513946e-16
13	-2.884235	2.15e-05	-0.0065599	7e-07	-0.0000930	-0.0023508	-0.0004230	0.0e+00	NA	-1e-07	0.0003408	NA	11	-3442.528	6907.101	71.98280	2.336900e-16
21	-2.881649	1.31e-05	-0.0065788	6e-07	0.0015150	-0.0013253	-0.0009513	1.8e-06	NA	-1e-07	NA	-0.0002092	11	-3442.806	6907.657	72.53841	1.770074e-16
7	-2.901278	1.25e-05	-0.0058625	8e-07	0.0016024	-0.0024539	0.0031392	1.8e-06	-0.0001449	-1e-07	NA	NA	11	-3443.046	6908.138	73.01974	1.391463e-16
30	-2.903249	2.15e-05	-0.0058572	7e-07	-0.0005994	-0.0010691	0.0037845	NA	-0.0001586	-1e-07	0.0004443	-0.0002653	12	-3444.762	6913.578	78.46031	9.163629e-18
29	-2.882138	1.85e-05	-0.0065971	6e-07	-0.0003202	-0.0010362	-0.0007469	7.0e-07	NA	-1e-07	0.0003770	-0.0002574	12	-3449.269	6922.592	87.47390	1.011093e-19
15	-2.903204	1.79e-05	-0.0058308	7e-07	-0.0000759	-0.0024175	0.0036222	7.0e-07	-0.0001538	-1e-07	0.0003472	NA	12	-3449.583	6923.220	88.10200	7.385839e-20
23	-2.905869	7.40e-06	-0.0056354	7e-07	0.0015654	-0.0011030	0.0042394	2.9e-06	-0.0002007	-1e-07	NA	-0.0002702	12	-3449.752	6923.557	88.43915	6.240085e-20
31	-2.909040	1.25e-05	-0.0055511	7e-07	-0.0003574	-0.0007748	0.0050239	1.9e-06	-0.0002227	-1e-07	0.0003961	-0.0003275	13	-3456.166	6938.395	103.27731	3.742176e-23

kable(dredged.ND.model, "html", caption = "ND Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

ND Dredge Table

	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	SDi	bioclim4:SDi	log(range.size.m2):SDi	NPP:SDi	PC1.LIG:SDi	residuals.PC1:SDi	df	logLik	AICc	delta	weight
0	0.1404205	5e-07	-0.0001462	0	-0.0000102	-6.47e-05	-0.0000575	NA	NA	NA	NA	NA	8	14329.64	-28643.26	0.00000	9.999550e-01
8	0.1403564	4e-07	-0.0001422	0	-0.0000850	-5.71e-05	-0.0000618	NA	NA	NA	1.48e-05	NA	9	14320.05	-28622.06	21.19677	2.495513e-05
16	0.1404275	5e-07	-0.0001466	0	-0.0000108	-5.67e-05	-0.0000581	NA	NA	NA	NA	-1.6e-06	9	14319.25	-28620.46	22.79683	1.121276e-05
2	0.1400429	5e-07	-0.0001322	0	-0.0000107	-6.61e-05	0.0000205	NA	-2.9e-06	NA	NA	NA	9	14318.99	-28619.95	23.30593	8.692847e-06
1	0.1403637	6e-07	-0.0001457	0	-0.0000101	-6.47e-05	-0.0000496	0e+00	NA	NA	NA	NA	9	14314.58	-28611.13	32.12698	1.056075e-07
4	0.1406803	5e-07	-0.0001499	0	-0.0000090	-6.52e-05	-0.0000932	NA	NA	0	NA	NA	9	14311.61	-28605.19	38.06105	5.434109e-09
24	0.1403735	4e-07	-0.0001429	0	-0.0000910	-3.30e-05	-0.0000638	NA	NA	NA	1.57e-05	-4.8e-06	10	14309.74	-28599.43	43.82308	3.047333e-10
10	0.1398108	4e-07	-0.0001218	0	-0.0000886	-5.88e-05	0.0000503	NA	-4.2e-06	NA	1.54e-05	NA	10	14309.46	-28598.88	44.37966	2.307068e-10
18	0.1400302	5e-07	-0.0001318	0	-0.0000116	-5.60e-05	0.0000242	NA	-3.1e-06	NA	NA	-2.1e-06	10	14308.61	-28597.18	46.08010	9.858583e-11
9	0.1401845	8e-07	-0.0001400	0	-0.0000977	-5.59e-05	-0.0000401	-1e-07	NA	NA	1.74e-05	NA	10	14305.25	-28590.46	52.79550	3.432282e-12
17	0.1403708	6e-07	-0.0001459	0	-0.0000104	-6.13e-05	-0.0000504	0e+00	NA	NA	NA	-7.0e-07	10	14304.23	-28588.42	54.83371	1.238774e-12
3	0.1400821	5e-07	-0.0001346	0	-0.0000106	-6.59e-05	0.0000102	0e+00	-2.4e-06	NA	NA	NA	10	14304.00	-28587.96	55.29503	9.835981e-13
12	0.1406450	4e-07	-0.0001461	0	-0.0000867	-5.74e-05	-0.0001020	NA	NA	0	1.54e-05	NA	10	14302.09	-28584.14	59.11661	1.455366e-13
20	0.1406914	5e-07	-0.0001503	0	-0.0000097	-5.54e-05	-0.0000943	NA	NA	0	NA	-2.0e-06	10	14301.22	-28582.41	60.84729	6.125736e-14
6	0.1404437	5e-07	-0.0001414	0	-0.0000094	-6.61e-05	-0.0000452	NA	-1.7e-06	0	NA	NA	10	14300.95	-28581.85	61.40235	4.641190e-14
26	0.1397605	4e-07	-0.0001201	0	-0.0000961	-3.08e-05	0.0000624	NA	-4.8e-06	NA	1.64e-05	-5.7e-06	11	14299.18	-28576.32	66.93341	2.921291e-15
5	0.1406261	6e-07	-0.0001493	0	-0.0000090	-6.53e-05	-0.0000857	0e+00	NA	0	NA	NA	10	14296.55	-28573.06	70.19982	5.705361e-16
25	0.1402082	8e-07	-0.0001405	0	-0.0000999	-4.30e-05	-0.0000430	-1e-07	NA	NA	1.76e-05	-2.6e-06	11	14294.92	-28567.80	75.45446	4.123369e-17
11	0.1399495	7e-07	-0.0001307	0	-0.0000978	-5.69e-05	0.0000099	-1e-07	-2.0e-06	NA	1.73e-05	NA	11	14294.66	-28567.27	75.98128	3.168506e-17
19	0.1400559	5e-07	-0.0001333	0	-0.0000113	-5.77e-05	0.0000174	0e+00	-2.7e-06	NA	NA	-1.7e-06	11	14293.70	-28565.35	77.90868	1.208721e-17
28	0.1406725	4e-07	-0.0001470	0	-0.0000935	-3.05e-05	-0.0001055	NA	NA	0	1.64e-05	-5.4e-06	11	14291.80	-28561.55	81.70210	1.813827e-18
14	0.1402430	4e-07	-0.0001316	0	-0.0000890	-5.86e-05	-0.0000209	NA	-3.0e-06	0	1.58e-05	NA	11	14291.46	-28560.87	82.38244	1.290804e-18
22	0.1404322	5e-07	-0.0001410	0	-0.0000103	-5.50e-05	-0.0000416	NA	-1.9e-06	0	NA	-2.3e-06	11	14290.57	-28559.08	84.17121	5.277578e-19
13	0.1404724	8e-07	-0.0001438	0	-0.0000991	-5.62e-05	-0.0000802	-1e-07	NA	0	1.79e-05	NA	11	14287.28	-28552.52	90.73241	1.984664e-20
21	0.1406404	6e-07	-0.0001497	0	-0.0000094	-5.93e-05	-0.0000874	0e+00	NA	0	NA	-1.2e-06	11	14286.20	-28550.36	92.89788	6.721422e-21
7	0.1404979	6e-07	-0.0001445	0	-0.0000092	-6.58e-05	-0.0000590	0e+00	-1.0e-06	0	NA	NA	11	14285.96	-28549.88	93.37861	5.285327e-21
27	0.1398914	7e-07	-0.0001279	0	-0.0001008	-3.94e-05	0.0000252	-1e-07	-2.8e-06	NA	1.76e-05	-3.6e-06	12	14284.39	-28544.73	98.52951	4.023185e-22
30	0.1401987	4e-07	-0.0001299	0	-0.0000970	-2.90e-05	-0.0000096	NA	-3.5e-06	0	1.69e-05	-6.0e-06	12	14281.20	-28538.34	104.91315	1.653402e-23
29	0.1405076	8e-07	-0.0001446	0	-0.0001019	-4.01e-05	-0.0000847	-1e-07	NA	0	1.82e-05	-3.3e-06	12	14276.97	-28529.89	113.36490	2.416050e-25
15	0.1404182	8e-07	-0.0001418	0	-0.0000991	-5.65e-05	-0.0000689	-1e-07	-4.0e-07	0	1.79e-05	NA	12	14276.70	-28529.34	113.91919	1.831227e-25
23	0.1404716	5e-07	-0.0001432	0	-0.0000099	-5.75e-05	-0.0000518	0e+00	-1.4e-06	0	NA	-1.7e-06	12	14275.66	-28527.26	115.99255	6.494083e-26
31	0.1403601	7e-07	-0.0001390	0	-0.0001023	-3.86e-05	-0.0000535	-1e-07	-1.2e-06	0	1.82e-05	-3.7e-06	13	14266.43	-28506.79	136.46456	2.328506e-30

From the dredged model list we can see that the top model is one with no interactions. Run this model on the MCC and 100 trees. Noting that each model uses a unique set of values for DR/ND and a unique tree in the correlation structure.

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
MCC.DR.top <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.DR.top, 'data/MCC.DR.top.rds')

#Run model for ND
MCC.ND.top <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.ND.top, 'data/MCC.ND.top.rds')
#Run the 100 models for DR and ND using the best model:

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ES")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
DR.model.data <- lapply(DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
DR.pgls.models <- mcmapply(function(x,y) {
  gls(ES ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(0.985, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = DR.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(DR.pgls.models, "data/DR.pgls.models.rds")

#Now for Node Density:
ND.model.data <- lapply(nd.list, function(x) {
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
ND.model.data <- lapply(ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(0.9996, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = ND.model.data, y = pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(ND.pgls.models, "data/ND.pgls.models.rds")

From these models we can get a distribution of the estimates for the models alongside the 95 % CIs of the MCC model, this will show the variation beween trees, and the variation associated with the top MCC model.

DR.pgls.models <- readRDS("data/DR.pgls.models.rds")
MCC.DR.top <- readRDS('data/MCC.DR.top.rds')

DR.pgls.summary <- lapply(DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.DR.summary <- data.frame(MCC.DR.top$coefficients, confint(MCC.DR.top)) %>% tibble::rownames_to_column()

DR.pgls.summary <- bind_rows(DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

DR.plot <-DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

DR.plot <- DR.plot + geom_errorbarh(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("DR Models")




#Now plot this for Node Density (ND):  

ND.pgls.models <- readRDS("data/ND.pgls.models.rds")
MCC.ND.top <- readRDS('data/MCC.ND.top.rds')

ND.pgls.summary <- lapply(ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.ND.summary <- data.frame(MCC.ND.top$coefficients, confint(MCC.ND.top)) %>% tibble::rownames_to_column()

ND.pgls.summary <- bind_rows(ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

ND.plot <-ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

ND.plot <- ND.plot + geom_errorbarh(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("ND Models")

PGLS Models on BAMM estimates

Using the previously generated BAMM dataframe we can now undertake model selection for speciation and extinction. The process is similar to that used on DR and ND, except given the use of a Bayesian approach in BAMM we can make use of varying levels of uncertainty between tips (species) by constructing a weighted model, where the weight is the inverse of the variance, meaning more precise estimates of speciation or extinction at a given tip (species) holds higher weight in the model.

Based on preliminary findings we found that Pagal’s lambda was = 1 and running corPagel lead to problemd of convergence. Therefore the following models are run assuming Brownian Motion with corBrownian.

MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')

#Create model dataframe for use in models
MCC.BAMM.model.data <- left_join(restricted.data %>% 
                                   dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
                                 MCC.BAMM.df %>% as.data.frame(), 
                                 by = "TipLabel")

saveRDS(MCC.BAMM.model.data, 'data/MCC.BAMM.model.data.rds')
#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.BAMM.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.BAMM.model.data) <- MCC.BAMM.model.data$TipLabel

#Run model for DR
MCC.BAMM.lambda <- gls(mean.lambda ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ sqrt(var.lambda), #sqrt to account for overdispersedskewed variance distribution
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

#Run model for DR
MCC.BAMM.mu <- gls(mean.mu ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")


#Dredge the global MCC models

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.BAMM.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.MCC.lambda <- pdredge(MCC.BAMM.lambda, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.lambda, "data/dredged.MCC.lambda.rds")

dredged.MCC.mu <- pdredge(MCC.BAMM.mu, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.mu, "data/dredged.MCC.mu.rds")

Table S6: The dredged models for BAMM speciation and BAMM extinction both show the top model is one with no interactions , with \(\delta AICc > 4\). This is the same situation as DR/ND models (see above).

dredged.MCC.lambda <- readRDS("data/dredged.MCC.lambda.rds")
dredged.MCC.mu <- readRDS("data/dredged.MCC.mu.rds")
kable(dredged.MCC.lambda, "html", caption = "BAMM-Speciation Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM-Speciation Dredge Table

	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	SDi	bioclim4:SDi	log(range.size.m2):SDi	NPP:SDi	PC1.LIG:SDi	residuals.PC1:SDi	df	logLik	AICc	delta	weight
0	-2.109226	-9.0e-07	-0.0002225	0	-0.0000487	0.0000043	-0.0000143	NA	NA	NA	NA	NA	8	6820.148	-13624.27	0.00000	9.998911e-01
16	-2.109176	-9.0e-07	-0.0002237	0	-0.0000508	0.0000683	-0.0000199	NA	NA	NA	NA	-1.38e-05	9	6811.014	-13604.00	20.27447	3.957376e-05
8	-2.109265	-1.0e-06	-0.0002207	0	-0.0000969	0.0000114	-0.0000161	NA	NA	NA	1.04e-05	NA	9	6810.971	-13603.91	20.35993	3.791842e-05
2	-2.109396	-9.0e-07	-0.0002160	0	-0.0000493	0.0000040	0.0000283	NA	-1.6e-06	NA	NA	NA	9	6810.770	-13603.51	20.76242	3.100641e-05
1	-2.109089	-1.2e-06	-0.0002250	0	-0.0000511	0.0000046	-0.0000299	1e-07	NA	NA	NA	NA	9	6806.364	-13594.70	29.57313	3.786430e-07
4	-2.108967	-9.0e-07	-0.0002264	0	-0.0000479	0.0000082	-0.0000514	NA	NA	0	NA	NA	9	6803.214	-13588.40	35.87346	1.622295e-08
24	-2.109217	-1.0e-06	-0.0002215	0	-0.0001185	0.0000950	-0.0000239	NA	NA	NA	1.46e-05	-1.74e-05	10	6801.896	-13583.75	40.51732	1.591214e-09
18	-2.109482	-9.0e-07	-0.0002120	0	-0.0000517	0.0000695	0.0000568	NA	-2.9e-06	NA	NA	-1.41e-05	10	6801.644	-13583.25	41.02023	1.237434e-09
10	-2.109495	-1.0e-06	-0.0002119	0	-0.0000984	0.0000111	0.0000414	NA	-2.2e-06	NA	1.06e-05	NA	10	6801.596	-13583.15	41.11709	1.178937e-09
17	-2.108913	-1.4e-06	-0.0002286	0	-0.0000555	0.0000827	-0.0000498	1e-07	NA	NA	NA	-1.68e-05	10	6797.294	-13574.55	49.72130	1.596283e-11
9	-2.109192	-1.1e-06	-0.0002221	0	-0.0000944	0.0000110	-0.0000239	0e+00	NA	NA	9.60e-06	NA	10	6797.223	-13574.41	49.86312	1.487013e-11
3	-2.109506	-1.3e-06	-0.0002080	0	-0.0000535	0.0000038	0.0000821	1e-07	-4.5e-06	NA	NA	NA	10	6797.094	-13574.15	50.12135	1.306895e-11
20	-2.108936	-9.0e-07	-0.0002273	0	-0.0000499	0.0000705	-0.0000544	NA	NA	0	NA	-1.34e-05	10	6794.078	-13568.12	56.15218	6.407108e-13
12	-2.108992	-1.0e-06	-0.0002247	0	-0.0000977	0.0000158	-0.0000554	NA	NA	0	1.08e-05	NA	10	6794.041	-13568.04	56.22755	6.170159e-13
6	-2.108998	-9.0e-07	-0.0002253	0	-0.0000480	0.0000081	-0.0000438	NA	-3.0e-07	0	NA	NA	10	6793.855	-13567.67	56.59936	5.123386e-13
26	-2.109645	-1.0e-06	-0.0002051	0	-0.0001220	0.0000975	0.0000831	NA	-4.0e-06	NA	1.50e-05	-1.80e-05	11	6792.534	-13563.02	61.24819	5.012765e-14
5	-2.108812	-1.2e-06	-0.0002292	0	-0.0000504	0.0000086	-0.0000693	1e-07	NA	0	NA	NA	10	6789.433	-13558.83	65.44243	6.156164e-15
25	-2.109027	-1.3e-06	-0.0002252	0	-0.0001141	0.0001021	-0.0000445	1e-07	NA	NA	1.29e-05	-1.91e-05	11	6788.188	-13554.33	69.94123	6.492436e-16
19	-2.109734	-1.5e-06	-0.0001945	0	-0.0000611	0.0000941	0.0001749	2e-07	-9.1e-06	NA	NA	-1.95e-05	11	6788.081	-13554.12	70.15475	5.835015e-16
11	-2.109552	-1.2e-06	-0.0002074	0	-0.0000951	0.0000101	0.0000733	1e-07	-3.9e-06	NA	9.30e-06	NA	11	6787.953	-13553.86	70.41057	5.134414e-16
28	-2.108963	-1.0e-06	-0.0002253	0	-0.0001189	0.0000977	-0.0000606	NA	NA	0	1.48e-05	-1.71e-05	11	6784.963	-13547.88	76.38995	2.582768e-17
22	-2.109125	-9.0e-07	-0.0002204	0	-0.0000506	0.0000710	-0.0000083	NA	-1.7e-06	0	NA	-1.37e-05	11	6784.727	-13547.41	76.86170	2.040073e-17
14	-2.109083	-1.0e-06	-0.0002214	0	-0.0000982	0.0000156	-0.0000334	NA	-8.0e-07	0	1.09e-05	NA	11	6784.683	-13547.32	76.95041	1.951567e-17
21	-2.108655	-1.4e-06	-0.0002325	0	-0.0000548	0.0000854	-0.0000867	1e-07	NA	0	NA	-1.65e-05	11	6780.361	-13538.68	85.59503	2.589570e-19
13	-2.108909	-1.2e-06	-0.0002264	0	-0.0000949	0.0000153	-0.0000645	0e+00	NA	0	9.90e-06	NA	11	6780.294	-13538.54	85.72902	2.421769e-19
7	-2.109119	-1.3e-06	-0.0002173	0	-0.0000521	0.0000078	0.0000103	1e-07	-3.0e-06	0	NA	NA	11	6780.178	-13538.31	85.95984	2.157800e-19
27	-2.109821	-1.5e-06	-0.0001922	0	-0.0001185	0.0001127	0.0001733	1e-07	-8.8e-06	NA	1.27e-05	-2.17e-05	12	6778.973	-13533.89	90.37795	2.369350e-20
30	-2.109279	-1.0e-06	-0.0002136	0	-0.0001213	0.0000993	0.0000163	NA	-2.8e-06	0	1.51e-05	-1.76e-05	12	6775.618	-13527.18	97.08814	8.270481e-22
29	-2.108760	-1.3e-06	-0.0002292	0	-0.0001143	0.0001051	-0.0000828	1e-07	NA	0	1.31e-05	-1.88e-05	12	6771.256	-13518.46	105.81172	1.054947e-23
23	-2.109417	-1.5e-06	-0.0002023	0	-0.0000598	0.0000946	0.0001146	2e-07	-7.8e-06	0	NA	-1.90e-05	12	6771.163	-13518.27	105.99910	9.605952e-24
15	-2.109145	-1.2e-06	-0.0002171	0	-0.0000953	0.0000146	-0.0000028	0e+00	-2.4e-06	0	9.70e-06	NA	12	6771.041	-13518.03	106.24362	8.500498e-24
31	-2.109486	-1.5e-06	-0.0002005	0	-0.0001180	0.0001135	0.0001092	1e-07	-7.4e-06	0	1.29e-05	-2.11e-05	13	6762.057	-13498.05	126.21895	3.907117e-28

kable(dredged.MCC.mu, "html", caption = "BAMM-Extinction Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM-Extinction Dredge Table

	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	SDi	bioclim4:SDi	log(range.size.m2):SDi	NPP:SDi	PC1.LIG:SDi	residuals.PC1:SDi	df	logLik	AICc	delta	weight
0	-4.355741	-1.3e-06	-0.0004220	0e+00	-0.0002276	-0.0000745	0.0000239	NA	NA	NA	NA	NA	8	2587.288	-5158.551	0.00000	9.997936e-01
16	-4.355884	-1.2e-06	-0.0004231	0e+00	-0.0002369	0.0000641	0.0000107	NA	NA	NA	NA	-2.98e-05	9	2578.877	-5139.723	18.82838	8.154169e-05
8	-4.355680	-1.3e-06	-0.0004247	0e+00	-0.0001838	-0.0000819	0.0000266	NA	NA	NA	-8.50e-06	NA	9	2578.674	-5139.317	19.23423	6.656558e-05
2	-4.356963	-1.3e-06	-0.0003779	0e+00	-0.0002336	-0.0000760	0.0002793	NA	-9.60e-06	NA	NA	NA	9	2578.528	-5139.025	19.52584	5.753445e-05
1	-4.355554	-1.7e-06	-0.0004251	0e+00	-0.0002285	-0.0000749	0.0000004	1e-07	NA	NA	NA	NA	9	2574.173	-5130.314	28.23710	7.384157e-07
4	-4.354974	-1.3e-06	-0.0004343	-1e-07	-0.0002234	-0.0000731	-0.0000587	NA	NA	0	NA	NA	9	2570.982	-5123.933	34.61801	3.038819e-08
24	-4.355849	-1.2e-06	-0.0004245	0e+00	-0.0002139	0.0000564	0.0000125	NA	NA	NA	-4.40e-06	-2.90e-05	10	2570.271	-5120.505	38.04631	5.473417e-09
18	-4.357443	-1.1e-06	-0.0003672	0e+00	-0.0002451	0.0000706	0.0003338	NA	-1.21e-05	NA	NA	-3.16e-05	10	2570.130	-5120.223	38.32863	4.752851e-09
10	-4.356868	-1.2e-06	-0.0003818	0e+00	-0.0001920	-0.0000829	0.0002739	NA	-9.30e-06	NA	-8.10e-06	NA	10	2569.914	-5119.791	38.76025	3.830269e-09
17	-4.355503	-2.0e-06	-0.0004301	0e+00	-0.0002406	0.0000909	-0.0000433	2e-07	NA	NA	NA	-3.58e-05	10	2565.829	-5111.620	46.93151	6.439636e-11
9	-4.355419	-1.8e-06	-0.0004295	0e+00	-0.0001730	-0.0000844	-0.0000034	1e-07	NA	NA	-1.08e-05	NA	10	2565.588	-5111.139	47.41245	5.063225e-11
3	-4.357282	-2.0e-06	-0.0003593	0e+00	-0.0002385	-0.0000775	0.0003758	2e-07	-1.49e-05	NA	NA	NA	10	2565.518	-5110.999	47.55279	4.720112e-11
20	-4.355085	-1.1e-06	-0.0004360	-1e-07	-0.0002327	0.0000686	-0.0000759	NA	NA	0	NA	-3.05e-05	10	2562.574	-5105.111	53.44028	2.485999e-12
12	-4.354937	-1.2e-06	-0.0004364	-1e-07	-0.0001843	-0.0000798	-0.0000544	NA	NA	0	-7.60e-06	NA	10	2562.368	-5104.698	53.85291	2.022555e-12
6	-4.355914	-1.2e-06	-0.0004017	-1e-07	-0.0002281	-0.0000743	0.0001311	NA	-6.80e-06	0	NA	NA	10	2562.236	-5104.434	54.11718	1.772201e-12
26	-4.357390	-1.1e-06	-0.0003691	0e+00	-0.0002265	0.0000643	0.0003302	NA	-1.20e-05	NA	-3.60e-06	-3.09e-05	11	2561.526	-5101.007	57.54451	3.193580e-13
5	-4.354737	-1.7e-06	-0.0004382	-1e-07	-0.0002242	-0.0000736	-0.0000879	1e-07	NA	0	NA	NA	10	2557.870	-5095.703	62.84860	2.251695e-14
19	-4.358289	-2.6e-06	-0.0003247	0e+00	-0.0002597	0.0001253	0.0005543	3e-07	-2.41e-05	NA	NA	-4.41e-05	11	2557.246	-5092.447	66.10399	4.421925e-15
25	-4.355407	-2.1e-06	-0.0004331	0e+00	-0.0002004	0.0000796	-0.0000449	2e-07	NA	NA	-7.80e-06	-3.48e-05	11	2557.245	-5092.444	66.10718	4.414887e-15
11	-4.357205	-2.1e-06	-0.0003614	0e+00	-0.0001787	-0.0000879	0.0003868	2e-07	-1.55e-05	NA	-1.18e-05	NA	11	2556.938	-5091.830	66.72172	3.246924e-15
28	-4.355066	-1.1e-06	-0.0004369	-1e-07	-0.0002154	0.0000628	-0.0000737	NA	NA	0	-3.30e-06	-2.98e-05	11	2553.970	-5085.894	72.65695	1.669761e-16
22	-4.356379	-1.1e-06	-0.0003913	-1e-07	-0.0002395	0.0000731	0.0001834	NA	-9.40e-06	0	NA	-3.18e-05	11	2553.839	-5085.633	72.91876	1.464880e-16
14	-4.355849	-1.2e-06	-0.0004048	-1e-07	-0.0001901	-0.0000806	0.0001292	NA	-6.60e-06	0	-7.40e-06	NA	11	2553.622	-5085.199	73.35222	1.179447e-16
21	-4.354615	-2.1e-06	-0.0004445	-1e-07	-0.0002364	0.0000981	-0.0001408	2e-07	NA	0	NA	-3.70e-05	11	2549.534	-5077.022	81.52954	1.976951e-18
13	-4.354629	-1.8e-06	-0.0004420	-1e-07	-0.0001726	-0.0000824	-0.0000896	1e-07	NA	0	-1.01e-05	NA	11	2549.286	-5076.526	82.02573	1.542585e-18
7	-4.356245	-2.0e-06	-0.0003832	-1e-07	-0.0002330	-0.0000758	0.0002278	2e-07	-1.21e-05	0	NA	NA	11	2549.225	-5076.405	82.14605	1.452519e-18
27	-4.358212	-2.7e-06	-0.0003270	0e+00	-0.0002158	0.0001133	0.0005585	4e-07	-2.43e-05	NA	-8.50e-06	-4.31e-05	12	2548.664	-5073.274	85.27712	3.035414e-19
30	-4.356346	-1.1e-06	-0.0003927	-1e-07	-0.0002249	0.0000682	0.0001817	NA	-9.30e-06	0	-2.80e-06	-3.12e-05	12	2545.236	-5066.417	92.13388	9.846754e-21
23	-4.357251	-2.6e-06	-0.0003486	-1e-07	-0.0002541	0.0001270	0.0004062	3e-07	-2.13e-05	0	NA	-4.41e-05	12	2540.954	-5057.854	100.69765	1.360485e-22
29	-4.354544	-2.1e-06	-0.0004469	-1e-07	-0.0002010	0.0000881	-0.0001407	2e-07	NA	0	-6.90e-06	-3.61e-05	12	2540.950	-5057.845	100.70611	1.354739e-22
15	-4.356211	-2.1e-06	-0.0003843	-1e-07	-0.0001773	-0.0000856	0.0002434	2e-07	-1.27e-05	0	-1.10e-05	NA	12	2540.644	-5057.235	101.31622	9.985545e-23
31	-4.357207	-2.6e-06	-0.0003500	-1e-07	-0.0002146	0.0001161	0.0004136	4e-07	-2.16e-05	0	-7.70e-06	-4.32e-05	13	2532.372	-5038.680	119.87105	9.337735e-27

Given that we are running models weighted according to the variance of the response (\(\lambda_{BAMM}\) and \(\mu_{BAMM}\)), we can check to see whether a weighted model isw favourable to an unweighted model using an anova to compare AIC values.

#Run the top model for the MCC
#Run model for ND
MCC.Lambda.top <- gls(mean.lambda ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top, 'data/MCC.Lambda.top.rds')

MCC.Mu.top <- gls(mean.mu ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

saveRDS(MCC.Mu.top, 'data/MCC.Mu.top.rds')


#We can also see how the models look without the weightings: 
MCC.Lambda.top.unweighted <- gls(mean.lambda ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top.unweighted, 'data/MCC.Lambda.top.unweighted.rds')


MCC.Mu.top.unweighted <- gls(mean.mu ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Mu.top.unweighted, 'data/MCC.Mu.top.unweighted.rds')

Table S6: Models with a set of weights for the response based on the inverse of the variance of the posterior distribution in \(\lambda_{BAMM}\) and \(\mu_{BAMM}\) have much lower AIC values, indicating that the weighting scheme improves the fit of the model. In any case, we found that an unweighted model did not change the main conclusions drawn in that they did not generate any statistically significant relationships.

MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')
MCC.Lambda.top.unweighted <- readRDS('data/MCC.Lambda.top.unweighted.rds')
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')
MCC.Mu.top.unweighted <- readRDS('data/MCC.Mu.top.unweighted.rds')

anova(MCC.Lambda.top, MCC.Lambda.top.unweighted) %>% pander(split.table = Inf)

	call	Model	df	AIC	BIC	logLik
MCC.Lambda.top	gls(model = mean.lambda ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~sqrt(var.lambda), method = “REML”)	1	8	-13624	-13571	6820
MCC.Lambda.top.unweighted	gls(model = mean.lambda ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”)	2	8	-11788	-11735	5902

anova(MCC.Mu.top, MCC.Mu.top.unweighted) %>% pander(split.table = Inf)

	call	Model	df	AIC	BIC	logLik
MCC.Mu.top	gls(model = mean.mu ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~sqrt(var.mu), method = “REML”)	1	8	-5159	-5105	2587
MCC.Mu.top.unweighted	gls(model = mean.mu ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”)	2	8	-5320	-5266	2668

We can run the top (no interactions) model on the 100 trees, each with unique estimates of speciation and extinction from the BAMM runs.

#Read in the BAMM data for the 100 trees
BAMM.df <- readRDS('data/BAMM.df.rds')

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the BAMM for the 100 trees

BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
BAMM.model.data <- lapply(BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.lambda.pgls.models, "data/BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.mu.pgls.models, "data/BAMM.mu.pgls.models.rds")

BAMM.lambda.pgls.models <- readRDS("data/BAMM.lambda.pgls.models.rds")
MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')

BAMM.lambda.pgls.summary <- lapply(BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

MCC.lambda.summary <- data.frame(MCC.Lambda.top$coefficients, confint(MCC.Lambda.top)) %>% tibble::rownames_to_column()

BAMM.lambda.pgls.summary <- bind_rows(BAMM.lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

  # for (x in unique(BAMM.lambda.pgls.summary$Parameter)[2:7]){
  #   filter(Parameter == x & between(Estimate, left = as.numeric(hpd.Lambda.top[1,x]), right = as.numeric(hpd.Lambda.top[2,x])))
  #   } 


remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  y
}

BAMM.lambda.pgls.summary.RO <- dcast(BAMM.lambda.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter)

## Using UCI as value column: use value.var to override.

BAMM.lambda.pgls.summary.RO$Estimate <- NULL
BAMM.lambda.pgls.summary.RO <- sapply(BAMM.lambda.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
BAMM.lambda.pgls.summary.RO <-melt(BAMM.lambda.pgls.summary.RO) %>% na.omit()
BAMM.lambda.pgls.summary.RO$Var1 <- NULL
colnames(BAMM.lambda.pgls.summary.RO) <- c("Parameter", "Estimate")

BAMM.lambda.plot <- BAMM.lambda.pgls.summary.RO %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.lambda.plot <- BAMM.lambda.plot + geom_errorbarh(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("BAMM Speciation")


#Figure for extinction

BAMM.mu.pgls.models <- readRDS("data/BAMM.mu.pgls.models.rds")
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')

BAMM.mu.pgls.summary <- lapply(BAMM.mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.mu.summary <- data.frame(MCC.Mu.top$coefficients, confint(MCC.Mu.top)) %>% tibble::rownames_to_column()

BAMM.mu.pgls.summary <- bind_rows(BAMM.mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

#Get rid of outliers
BAMM.mu.pgls.summary.RO <- dcast(BAMM.mu.pgls.summary %>% select(Parameter,Estimate) %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter)

## Using Estimate as value column: use value.var to override.

BAMM.mu.pgls.summary.RO$Estimate <- NULL
BAMM.mu.pgls.summary.RO <- sapply(BAMM.mu.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
BAMM.mu.pgls.summary.RO <-melt(BAMM.mu.pgls.summary.RO) %>% na.omit()
BAMM.mu.pgls.summary.RO$Var1 <- NULL
colnames(BAMM.mu.pgls.summary.RO) <- c("Parameter", "Estimate")

BAMM.mu.plot <-BAMM.mu.pgls.summary.RO %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.mu.plot <- BAMM.mu.plot + geom_errorbarh(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,0.5,.5),"cm"))+
  ggtitle("BAMM Extinction")

PGLS Model Summary

grid.arrange(
  symmetrise_scale(DR.plot, "x"),
  symmetrise_scale(ND.plot, "x"),
  symmetrise_scale(BAMM.lambda.plot, "x"),
  symmetrise_scale(BAMM.mu.plot, "x"),
  nrow = 1
)

Figure S8: This figure is the same basic figure as seen in the manuscript (FIgure 1), it provides model estimates for four response variables across 100 random trees alongside the MCC tree.

Table S7: The estimates of the MCC models plotted above are based on the following data tables.

MCC.DR.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates")

MCC DR Estimates

Parameter	Estimate	LCI	UCI
SDi	-0.001279	-0.003324	0.0007665
log(range.size.m2)	-0.006578	-0.01078	-0.002381
bioclim4	2.141e-05	-3.369e-05	7.65e-05
residuals.PC1	-0.00243	-0.007562	0.002703
PC1.LIG	0.001586	-0.003676	0.006849
NPP	2.937e-07	-1.321e-06	1.909e-06

MCC.ND.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates")

MCC ND Estimates

Parameter	Estimate	LCI	UCI
SDi	-5.745e-05	-0.0001239	8.957e-06
log(range.size.m2)	-0.0001462	-0.0002827	-9.773e-06
bioclim4	4.523e-07	-1.426e-06	2.33e-06
residuals.PC1	-6.469e-05	-0.0002368	0.0001074
PC1.LIG	-1.017e-05	-0.0001763	0.000156
NPP	-9.069e-09	-6.454e-08	4.64e-08

MCC.lambda.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates")

MCC BAMM Speciation Estimates

Parameter	Estimate	LCI	UCI
SDi	-1.43e-05	-0.0002574	0.0002288
log(range.size.m2)	-0.0002225	-0.00069	0.000245
bioclim4	-9.418e-07	-7.242e-06	5.359e-06
residuals.PC1	4.311e-06	-0.0005798	0.0005884
PC1.LIG	-4.874e-05	-0.0005747	0.0004772
NPP	-2.09e-08	-2.282e-07	1.864e-07

MCC.mu.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates")

MCC BAMM Extinction Estimates

Parameter	Estimate	LCI	UCI
SDi	2.39e-05	-0.0004098	0.0004576
log(range.size.m2)	-0.000422	-0.00133	0.0004858
bioclim4	-1.289e-06	-1.401e-05	1.143e-05
residuals.PC1	-7.453e-05	-0.001231	0.001082
PC1.LIG	-0.0002276	-0.001284	0.0008287
NPP	-2.024e-08	-4.178e-07	3.773e-07

From the estimates of 100 trees we can calculate HPD intervals. Note that these do no account for the CI of each estimate. Thus it is a more a representation of tree uncertainty than model uncertainty.

Table S8: Highest posterior density intervals are calculated for the above figure (ie; models using RGB values of sexual dichromatism). The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -).

hpd.DR.top <- list()
for(x in unique(DR.pgls.summary$Parameter)) {
hpd.DR.top[[x]] = hdi(DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select("Estimate"))
}
hpd.DR.top <- bind_rows(hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.DR.top, 'data/hpd.DR.top.rds')

#For ND
hpd.ND.top <- list()
for(x in unique(ND.pgls.summary$Parameter)) {
hpd.ND.top[[x]] = hdi(ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.ND.top <- bind_rows(hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 

saveRDS(hpd.ND.top, 'data/hpd.ND.top.rds')

hpd.Lambda.top <- list()
for(x in unique(BAMM.lambda.pgls.summary$Parameter)) {
hpd.Lambda.top[[x]] = hdi(BAMM.lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Lambda.top <- bind_rows(hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 


saveRDS(hpd.Lambda.top, 'data/hpd.Lambda.top.rds')

hpd.Mu.top <- list()
for(x in unique(BAMM.mu.pgls.summary$Parameter)) {
hpd.Mu.top[[x]] = hdi(BAMM.mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Mu.top <- bind_rows(hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.Mu.top, 'data/hpd.Mu.top.rds')

hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "DR HPD Interval")

DR HPD Interval

	SDi	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.00193	-0.00844	-2.16e-05	-0.00323	-0.00389	-1.4e-06
Upper	0.00156	-0.00182	5.6e-05	0.00479	0.00523	1.53e-06

hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "ND HPD Interval")

ND HPD Interval

	SDi	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-7.68e-05	-0.00019	-9.72e-07	-9.09e-05	-0.00015	-4.27e-08
Upper	5.97e-05	9.11e-06	1.62e-06	0.000175	0.000117	4.28e-08

hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Speciation HPD Interval")

BAMM Speciation HPD Interval

	SDi	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.00831	-0.00786	-0.000169	-0.0146	-0.0107	-5.56e-06
Upper	0.00474	0.0114	0.000124	0.0185	0.0084	2.42e-06

hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Extinction HPD Interval")

BAMM Extinction HPD Interval

	SDi	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.00788	-0.0165	-0.000159	-0.0216	-0.0161	-3.96e-06
Upper	0.0121	0.0237	0.000406	0.0367	0.0264	9.83e-06

From the HPD 95 % intervals and the intervals for the 95 % CI we can gain an estimate of how much variability there is within a model in comparison to variability between trees. To do this we calculate the width of the HPD interval relativel to the 95 % CI interval from the MCC.

Table S9: Model estimates from \(\lambda_{BAMM}\) and \(\mu_{BAMM}\) tip-rate estimates show high levels of variability between trees relative to estimates from \(\lambda_{DR}\) and \(\lambda_{ND}\). The interval ratio is calculated as the wisth of the HPD 95 % CI relative to the MCC 95 % CI. A value of 1 suggests that model estimates from 95 % of the trees fall within the MCC 95 % CI. Here, the large values for the estimates from BAMM tip-rates suggest that there is high variation in tip-rate estimates across trees that is not see within the DR and ND models.

hpd.DR.top2 <- hpd.DR.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.DR.top2$HPD.Interval <- hpd.DR.top2$Upper - hpd.DR.top2$Lower
DR.intervals <- left_join(hpd.DR.top2, MCC.DR.summary, by = "Parameter")
DR.intervals$MCC.Interval <- DR.intervals$UCI - DR.intervals$LCI
DR.intervals$DR.IntervalRatio <- hpd.DR.top2$HPD.Interval/DR.intervals$MCC.Interval

hpd.ND.top2 <- hpd.ND.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.ND.top2$HPD.Interval <- hpd.ND.top2$Upper - hpd.ND.top2$Lower
ND.intervals <- left_join(hpd.ND.top2, MCC.ND.summary, by = "Parameter")
ND.intervals$MCC.Interval <- ND.intervals$UCI - ND.intervals$LCI
ND.intervals$ND.IntervalRatio <- hpd.ND.top2$HPD.Interval/ND.intervals$MCC.Interval


hpd.Lambda.top2 <- hpd.Lambda.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.Lambda.top2$HPD.Interval <- hpd.Lambda.top2$Upper - hpd.Lambda.top2$Lower
Lambda.intervals <- left_join(hpd.Lambda.top2, MCC.lambda.summary, by = "Parameter")
Lambda.intervals$MCC.Interval <- Lambda.intervals$UCI - Lambda.intervals$LCI
Lambda.intervals$Lambda.IntervalRatio <- hpd.Lambda.top2$HPD.Interval/Lambda.intervals$MCC.Interval


hpd.Mu.top2 <- hpd.Mu.top %>% t() %>% as.data.frame() %>% rownames_to_column(var = "Parameter")
hpd.Mu.top2$HPD.Interval <- hpd.Mu.top2$Upper - hpd.Mu.top2$Lower
Mu.intervals <- left_join(hpd.Mu.top2, MCC.mu.summary, by = "Parameter")
Mu.intervals$MCC.Interval <- Mu.intervals$UCI - Mu.intervals$LCI
Mu.intervals$Mu.IntervalRatio <- hpd.Mu.top2$HPD.Interval/Mu.intervals$MCC.Interval


plyr::join_all(list(DR.intervals %>% select(Parameter, DR.IntervalRatio),
          ND.intervals %>% select(Parameter, ND.IntervalRatio),
          Lambda.intervals %>% select(Parameter, Lambda.IntervalRatio),
          Mu.intervals %>% select(Parameter, Mu.IntervalRatio)), by = "Parameter", type = "left") %>% `colnames<-`(c("Parameter", "λDR.Interval.Ratio", "λND.Interval.Ratio", "λBAMM.Interval.Ratio", "μBAMM.Interval.Ratio")) %>%
  pander(digits = 3, split.table = Inf) 

Parameter	λDR.Interval.Ratio	λND.Interval.Ratio	λBAMM.Interval.Ratio	μBAMM.Interval.Ratio
SDi	0.855	1.03	26.8	23.1
log(range.size.m2)	0.788	0.729	20.6	22.2
bioclim4	0.704	0.691	23.3	22.2
residuals.PC1	0.782	0.774	28.3	25.2
PC1.LIG	0.866	0.805	18.1	20.1
NPP	0.906	0.771	19.3	17.3

Subsetted analysis with Armenta data

Using the dataset from Armenta et al. (2008) we can conduct a subsetted analysis where dichromatism values are measured from spectrophotmetry and synthesised into a colour discriminability measure, which is a difference measure between the sexes. To make this dataset comprable with ours we use the absolute difference between the sexes. Therefore making the scale from monochromatism to dichromatism rather than female colouration to male colouration.

#Read in Armenta data
Armenta.data <- read.csv('data/Armenta_2008.csv', stringsAsFactors = F)
MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#A couple of rows have "-", remove from dataset
Armenta.data <- Armenta.data %>% dplyr::select(binomial, Colour.discriminability) %>% mutate(Colour.discriminability = replace(Colour.discriminability, Colour.discriminability == "-", "NA")) %>% filter(Colour.discriminability != "NA")
Armenta.data$Colour.discriminability <- Armenta.data$Colour.discriminability %>% as.numeric()

MCC.Armenta.model.data <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), Armenta.data, by = "binomial")

MCC.Armenta.model.data <- inner_join(MCC.Armenta.model.data, MCC.BAMM.df, by = "TipLabel") %>% filter(Colour.discriminability != "NA")

#Plot the correlation
p1 <- MCC.Armenta.model.data %>% ggplot(aes(x = SDi, y = Colour.discriminability))+
  geom_point()+
  geom_smooth(method = 'loess')+
  theme_minimal()

grid.newpage()
ggExtra::ggMarginal(
  p = p1,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

Figure S9: There is a correlation between the RGB measures and the Spectophotometry measures, The RGB measures seem to be more noisy around the lower values.

R-squared for relationship between RGB and Spectrophotometry Data:

summary(lm(SDi ~ Colour.discriminability,
   data = MCC.Armenta.model.data))$r.squared %>% pander()

0.6166

cor(MCC.Armenta.model.data$SDi, MCC.Armenta.model.data$Colour.discriminability)

## [1] 0.785264

Run PGLS model for the data

#Prune tree
pruned.MCC.Armenta.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.Armenta.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.Armenta.model.data) <- MCC.Armenta.model.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Armenta.global <- gls(MCC.DR ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Armenta.global <- gls(MCC.ND ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Armenta.global <- gls(mean.lambda ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Armenta.global <- gls(mean.mu ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Colour.discriminability*log(range.size.m2)
                         + Colour.discriminability*bioclim4
                         + Colour.discriminability*residuals.PC1
                         + Colour.discriminability*PC1.LIG
                         + Colour.discriminability*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.Armenta.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Armenta.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Armenta.dredged.ND.model <- pdredge(MCC.ND.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.ND.model, 'data/Armenta.dredged.ND.model.rds')

Armenta.dredged.DR.model <- pdredge(MCC.DR.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.DR.model, 'data/Armenta.dredged.DR.model.rds')

Armenta.dredged.spec.model <- pdredge(MCC.Lambda.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.spec.model, 'data/Armenta.dredged.spec.model.rds')

Armenta.dredged.extinct.model <- pdredge(MCC.Extinction.Armenta.global, fixed = c("Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.extinct.model, 'data/Armenta.dredged.extinct.model.rds')

Table S10: The dredged models for Models using a subset dataset where sexual dichromatism was measured using spectrophotometry, all four show the top model is one with no interactions , with \(\delta AICc > 4\).

Armenta.dredged.ND.model <- readRDS("data/Armenta.dredged.ND.model.rds")
Armenta.dredged.DR.model <- readRDS("data/Armenta.dredged.DR.model.rds")
Armenta.dredged.spec.model <- readRDS("data/Armenta.dredged.spec.model.rds")
Armenta.dredged.extinct.model <- readRDS("data/Armenta.dredged.extinct.model.rds")


kable(Armenta.dredged.ND.model, "html", caption = "ND Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

ND Spectrophotometry Dichromatism Dredge Table

	(Intercept)	bioclim4	Colour.discriminability	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	bioclim4:Colour.discriminability	Colour.discriminability:log(range.size.m2)	Colour.discriminability:NPP	Colour.discriminability:PC1.LIG	Colour.discriminability:residuals.PC1	df	logLik	AICc	delta	weight
0	0.2077881	6.9e-06	0.0000828	-0.0000378	5e-07	-0.0008841	0.0010532	NA	NA	NA	NA	NA	8	1033.2946	-2050.337	0.00000	9.972312e-01
8	0.2128993	6.2e-06	0.0003666	-0.0002044	4e-07	-0.0008428	0.0010460	NA	NA	NA	-0.0009029	NA	9	1027.8860	-2037.457	12.88059	1.591519e-03
2	0.2060699	7.3e-06	-0.0081861	0.0000182	5e-07	-0.0008841	0.0010718	NA	0.0002887	NA	NA	NA	9	1027.0481	-2035.781	14.55654	6.884668e-04
16	0.2083870	6.9e-06	0.0001261	-0.0000596	5e-07	-0.0008744	0.0010533	NA	NA	NA	NA	-0.0001539	9	1026.6890	-2035.063	15.27474	4.807590e-04
1	0.2072627	7.1e-06	-0.0002087	-0.0000199	5e-07	-0.0008882	0.0010535	6.00e-07	NA	NA	NA	NA	9	1021.9461	-2025.577	24.76044	4.189235e-06
10	0.2094193	7.0e-06	-0.0259111	-0.0000910	5e-07	-0.0008267	0.0011025	NA	0.0009215	NA	-0.0012568	NA	10	1022.3412	-2024.296	26.04102	2.208311e-06
24	0.2127267	6.2e-06	0.0003401	-0.0001949	4e-07	-0.0008524	0.0010447	NA	NA	NA	-0.0010503	0.0002588	10	1021.4594	-2022.533	27.80461	9.143286e-07
4	0.2077525	7.1e-06	0.0008792	-0.0000384	5e-07	-0.0009151	0.0010575	NA	NA	-1e-07	NA	NA	9	1019.5502	-2020.785	29.55223	3.816023e-07
18	0.2066684	7.2e-06	-0.0081477	-0.0000036	5e-07	-0.0008744	0.0010719	NA	0.0002889	NA	NA	-0.0001540	10	1020.4432	-2020.500	29.83700	3.309587e-07
9	0.2095915	7.9e-06	-0.0054207	-0.0000803	4e-07	-0.0008660	0.0010408	1.38e-05	NA	NA	-0.0023002	NA	10	1018.8039	-2017.222	33.11555	6.424583e-08
3	0.2062134	7.2e-06	-0.0088719	0.0000131	5e-07	-0.0008816	0.0010736	-4.00e-07	0.0003190	NA	NA	NA	10	1015.8364	-2011.287	39.05055	3.304277e-09
17	0.2074141	7.1e-06	-0.0005561	-0.0000272	5e-07	-0.0008795	0.0010540	1.60e-06	NA	NA	NA	-0.0002296	10	1015.4839	-2010.582	39.75568	2.322516e-09
26	0.2084783	7.2e-06	-0.0307006	-0.0000533	5e-07	-0.0008411	0.0011104	NA	0.0010878	NA	-0.0015881	0.0004694	11	1016.1516	-2009.839	40.49816	1.602260e-09
12	0.2133899	6.8e-06	0.0028641	-0.0002261	4e-07	-0.0009338	0.0010583	NA	NA	-3e-07	-0.0010090	NA	10	1014.3466	-2008.307	42.03013	7.448502e-10
6	0.2060643	7.3e-06	-0.0082262	0.0000185	5e-07	-0.0008837	0.0010718	NA	0.0002898	0e+00	NA	NA	10	1013.4370	-2006.488	43.84935	2.999366e-10
20	0.2083692	7.1e-06	0.0009782	-0.0000611	5e-07	-0.0009072	0.0010578	NA	NA	-1e-07	NA	-0.0001591	10	1012.9494	-2005.513	44.82470	1.841768e-10
11	0.2078948	8.2e-06	-0.0201410	-0.0000266	4e-07	-0.0008543	0.0010742	1.24e-05	0.0005364	NA	-0.0023673	NA	11	1012.8327	-2003.202	47.13593	5.799078e-11
25	0.2095645	7.9e-06	-0.0053564	-0.0000780	4e-07	-0.0008695	0.0010404	1.36e-05	NA	NA	-0.0023418	0.0001047	11	1012.3028	-2002.142	48.19581	3.413571e-11
19	0.2065718	7.2e-06	-0.0071738	-0.0000002	5e-07	-0.0008763	0.0010694	6.00e-07	0.0002463	NA	NA	-0.0001812	11	1009.3819	-1996.300	54.03762	1.839399e-12
5	0.2074923	7.2e-06	0.0006301	-0.0000293	5e-07	-0.0009132	0.0010571	3.00e-07	NA	-1e-07	NA	NA	10	1008.2734	-1996.161	54.17659	1.715926e-12
14	0.2096468	7.1e-06	-0.0243990	-0.0000994	4e-07	-0.0008428	0.0011022	NA	0.0008833	0e+00	-0.0012602	NA	11	1008.7339	-1995.004	55.33369	9.621384e-13
28	0.2132308	6.8e-06	0.0030279	-0.0002169	4e-07	-0.0009518	0.0010578	NA	NA	-3e-07	-0.0011857	0.0002957	11	1007.9537	-1993.443	56.89410	4.409596e-13
22	0.2067012	7.2e-06	-0.0079237	-0.0000048	5e-07	-0.0008768	0.0010718	NA	0.0002832	0e+00	NA	-0.0001544	11	1006.8344	-1991.205	59.13254	1.439887e-13
27	0.2074778	8.2e-06	-0.0230409	-0.0000101	4e-07	-0.0008604	0.0010801	1.17e-05	0.0006477	NA	-0.0024811	0.0002517	12	1006.4510	-1988.353	61.98466	3.459381e-14
13	0.2096921	7.9e-06	-0.0050367	-0.0000845	4e-07	-0.0008772	0.0010424	1.36e-05	NA	0e+00	-0.0022975	NA	11	1005.1167	-1987.770	62.56792	2.584311e-14
7	0.2062406	7.2e-06	-0.0087170	0.0000121	5e-07	-0.0008834	0.0010736	-4.00e-07	0.0003153	0e+00	NA	NA	11	1002.2404	-1982.017	68.32055	1.456052e-15
21	0.2075602	7.2e-06	0.0000142	-0.0000332	5e-07	-0.0008964	0.0010563	1.30e-06	NA	-1e-07	NA	-0.0002203	11	1001.8148	-1981.166	69.17175	9.513509e-16
30	0.2086525	7.3e-06	-0.0295459	-0.0000597	5e-07	-0.0008532	0.0011102	NA	0.0010585	0e+00	-0.0015897	0.0004681	12	1002.5432	-1980.537	69.80031	6.947876e-16
15	0.2073335	8.1e-06	-0.0234976	-0.0000057	4e-07	-0.0008176	0.0010741	1.27e-05	0.0006154	1e-07	-0.0023857	NA	12	999.2588	-1973.968	76.36914	2.602837e-17
29	0.2096993	7.9e-06	-0.0048284	-0.0000836	4e-07	-0.0008852	0.0010425	1.33e-05	NA	0e+00	-0.0023414	0.0001133	12	998.6289	-1972.708	77.62895	1.386383e-17
23	0.2065684	7.2e-06	-0.0071929	-0.0000001	5e-07	-0.0008760	0.0010694	6.00e-07	0.0002468	0e+00	NA	-0.0001813	12	995.7873	-1967.025	83.31204	8.087520e-19
31	0.2069328	8.1e-06	-0.0262977	0.0000102	4e-07	-0.0008246	0.0010799	1.20e-05	0.0007241	1e-07	-0.0024983	0.0002501	13	992.8768	-1959.112	91.22581	1.546546e-20

kable(Armenta.dredged.DR.model, "html", caption = "DR Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

DR Spectrophotometry Dichromatism Dredge Table

	(Intercept)	bioclim4	Colour.discriminability	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	bioclim4:Colour.discriminability	Colour.discriminability:log(range.size.m2)	Colour.discriminability:NPP	Colour.discriminability:PC1.LIG	Colour.discriminability:residuals.PC1	df	logLik	AICc	delta	weight
0	-2.677136	0.0001789	0.0070708	0.0090337	7.8e-06	-0.0315655	0.0156001	NA	NA	NA	NA	NA	9	-512.3837	1043.083	0.000000	9.479871e-01
8	-2.630812	0.0001604	0.0092806	0.0078425	7.2e-06	-0.0288585	0.0166944	NA	NA	NA	-0.0169326	NA	10	-514.7327	1049.851	6.768670	3.213693e-02
2	-2.669309	0.0001769	0.0667921	0.0087960	7.7e-06	-0.0314417	0.0155465	NA	-0.0020932	NA	NA	NA	10	-515.8406	1052.067	8.984518	1.061303e-02
16	-2.675941	0.0001777	0.0078042	0.0089905	7.8e-06	-0.0313098	0.0157537	NA	NA	NA	NA	-0.0017242	10	-516.0987	1052.583	9.500619	8.199156e-03
10	-2.653060	0.0001640	-0.2224726	0.0085957	7.3e-06	-0.0288185	0.0171615	NA	0.0081372	NA	-0.0199731	NA	11	-517.9245	1058.313	15.230367	4.672726e-04
24	-2.624359	0.0001602	0.0064975	0.0077898	6.9e-06	-0.0292878	0.0163567	NA	NA	NA	-0.0213744	0.0079966	11	-518.0589	1058.582	15.499006	4.085413e-04
18	-2.668225	0.0001757	0.0668268	0.0087564	7.8e-06	-0.0311887	0.0157002	NA	-0.0020689	NA	NA	-0.0017135	11	-519.5546	1061.573	18.490405	9.155078e-05
1	-2.671939	0.0001754	0.0136694	0.0088999	7.8e-06	-0.0312585	0.0157116	-0.0000153	NA	NA	NA	NA	10	-520.8980	1062.182	19.099199	6.752492e-05
9	-2.642496	0.0001826	-0.0714208	0.0081999	7.2e-06	-0.0295284	0.0166390	0.0001929	NA	NA	-0.0366791	NA	11	-521.5732	1065.610	22.527776	1.216068e-05
4	-2.673868	0.0001865	0.0487209	0.0088595	7.7e-06	-0.0331665	0.0160592	NA	NA	-4.6e-06	NA	NA	10	-523.1346	1066.655	23.572501	7.212717e-06
26	-2.654133	0.0001654	-0.3246079	0.0088415	7.0e-06	-0.0293710	0.0169056	NA	0.0115950	NA	-0.0270764	0.0104643	12	-521.0596	1066.669	23.585830	7.164808e-06
3	-2.669100	0.0001751	0.0435357	0.0088140	7.7e-06	-0.0312622	0.0156661	-0.0000116	-0.0011029	NA	NA	NA	11	-524.1891	1070.842	27.759554	8.889798e-07
17	-2.672419	0.0001753	0.0125561	0.0089088	7.8e-06	-0.0311560	0.0157961	-0.0000115	NA	NA	NA	-0.0011657	11	-524.5318	1071.528	28.444956	6.310427e-07
12	-2.618510	0.0001699	0.0766749	0.0073776	6.8e-06	-0.0310209	0.0175985	NA	NA	-7.4e-06	-0.0195276	NA	11	-525.1224	1072.709	29.626063	3.496106e-07
11	-2.643111	0.0001824	-0.0776207	0.0082299	7.2e-06	-0.0295154	0.0166688	0.0001919	0.0002323	NA	-0.0366674	NA	12	-524.8734	1074.296	31.213509	1.580793e-07
25	-2.636482	0.0001821	-0.0727080	0.0081501	6.9e-06	-0.0299139	0.0163324	0.0001898	NA	NA	-0.0404677	0.0073866	12	-524.9401	1074.429	31.346776	1.478892e-07
6	-2.638321	0.0001829	0.3290294	0.0077113	7.3e-06	-0.0335589	0.0160458	NA	-0.0090012	-7.2e-06	NA	NA	11	-526.2544	1074.973	31.890077	1.127094e-07
20	-2.672935	0.0001854	0.0488811	0.0088282	7.7e-06	-0.0329441	0.0161800	NA	NA	-4.6e-06	NA	-0.0013734	11	-526.8525	1076.169	33.086302	6.197310e-08
19	-2.668293	0.0001750	0.0571451	0.0087749	7.7e-06	-0.0311417	0.0157386	-0.0000048	-0.0016595	NA	NA	-0.0014888	12	-527.7809	1080.111	37.028433	8.633371e-09
28	-2.607544	0.0001714	0.0853530	0.0072208	6.4e-06	-0.0320005	0.0173071	NA	NA	-8.8e-06	-0.0259988	0.0107877	12	-528.2316	1081.012	37.929758	5.501233e-09
14	-2.622604	0.0001699	0.0367068	0.0075259	6.9e-06	-0.0309141	0.0176481	NA	0.0013029	-7.1e-06	-0.0199055	NA	12	-528.3498	1081.249	38.166275	4.887662e-09
27	-2.644834	0.0001822	-0.1679306	0.0084555	7.0e-06	-0.0298937	0.0165119	0.0001781	0.0035063	NA	-0.0410063	0.0081624	13	-528.1672	1082.976	39.893793	2.060513e-09
22	-2.637873	0.0001820	0.3270864	0.0076981	7.3e-06	-0.0333773	0.0161453	NA	-0.0089347	-7.2e-06	NA	-0.0010874	12	-529.9733	1084.496	41.413285	9.638776e-10
5	-2.656263	0.0001785	0.0857028	0.0083648	7.4e-06	-0.0328347	0.0165795	-0.0000478	NA	-6.4e-06	NA	NA	11	-531.4166	1085.297	42.214478	6.457213e-10
13	-2.635478	0.0001841	-0.0303738	0.0079455	7.0e-06	-0.0304341	0.0170601	0.0001683	NA	-3.4e-06	-0.0353444	NA	12	-532.3431	1089.235	46.152785	9.012648e-11
30	-2.620812	0.0001720	-0.0499593	0.0076922	6.5e-06	-0.0317146	0.0174172	NA	0.0044218	-7.8e-06	-0.0276394	0.0114010	13	-531.4062	1089.454	46.371621	8.078539e-11
7	-2.635314	0.0001785	0.2873035	0.0076666	7.2e-06	-0.0332479	0.0163926	-0.0000308	-0.0068943	-7.8e-06	NA	NA	12	-534.5375	1093.624	50.541586	1.004237e-11
21	-2.654840	0.0001782	0.0907471	0.0083515	7.3e-06	-0.0330081	0.0165591	-0.0000543	NA	-6.7e-06	NA	0.0014471	12	-535.0039	1094.557	51.474387	6.299140e-12
15	-2.627129	0.0001839	0.0527643	0.0076688	6.9e-06	-0.0306250	0.0169889	0.0001727	-0.0027968	-4.0e-06	-0.0349658	NA	13	-535.5408	1097.724	54.640887	1.293256e-12
29	-2.624974	0.0001842	-0.0138454	0.0077690	6.6e-06	-0.0313081	0.0168654	0.0001537	NA	-4.9e-06	-0.0393974	0.0090385	13	-535.5861	1097.814	54.731532	1.235951e-12
23	-2.635166	0.0001782	0.2848941	0.0076832	7.2e-06	-0.0333049	0.0164014	-0.0000341	-0.0067339	-7.9e-06	NA	0.0006111	13	-538.1141	1102.870	59.787439	9.865632e-14
31	-2.625387	0.0001840	-0.0182030	0.0077919	6.6e-06	-0.0312931	0.0168860	0.0001531	0.0001569	-4.9e-06	-0.0393997	0.0090593	14	-538.7716	1106.285	63.202639	1.788643e-14

kable(Armenta.dredged.spec.model, "html", caption = "BAMM Speciation Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM Speciation Spectrophotometry Dichromatism Dredge Table

	(Intercept)	bioclim4	Colour.discriminability	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	bioclim4:Colour.discriminability	Colour.discriminability:log(range.size.m2)	Colour.discriminability:NPP	Colour.discriminability:PC1.LIG	Colour.discriminability:residuals.PC1	df	logLik	AICc	delta	weight
0	-2.482091	-2.59e-05	-0.0005285	0.0018448	-1.6e-06	-0.0042436	0.0024890	NA	NA	NA	NA	NA	8	181.4767	-346.7017	0.00000	9.919991e-01
2	-2.488239	-2.42e-05	-0.0427556	0.0020142	-1.6e-06	-0.0042772	0.0025602	NA	0.0014676	NA	NA	NA	9	176.7474	-335.1796	11.52204	3.122707e-03
8	-2.477350	-2.71e-05	-0.0001200	0.0017133	-1.7e-06	-0.0041157	0.0025958	NA	NA	NA	-0.0011826	NA	9	176.5780	-334.8407	11.86100	2.635893e-03
16	-2.478790	-2.64e-05	-0.0000802	0.0017377	-1.6e-06	-0.0040974	0.0025735	NA	NA	NA	NA	-0.0009394	9	176.3953	-334.4753	12.22638	2.195774e-03
1	-2.475060	-2.79e-05	0.0040917	0.0016355	-1.7e-06	-0.0041107	0.0025852	-1.04e-05	NA	NA	NA	NA	9	171.7160	-325.1168	21.58489	2.038976e-05
10	-2.484288	-2.53e-05	-0.0750059	0.0019038	-1.6e-06	-0.0040674	0.0028130	NA	0.0026147	NA	-0.0021843	NA	10	172.1264	-323.8668	22.83488	1.091389e-05
18	-2.484933	-2.47e-05	-0.0422043	0.0019071	-1.6e-06	-0.0041314	0.0026442	NA	0.0014640	NA	NA	-0.0009365	10	171.6660	-322.9461	23.75558	6.887368e-06
24	-2.476734	-2.71e-05	0.0000344	0.0016866	-1.7e-06	-0.0040638	0.0026159	NA	NA	NA	-0.0008273	-0.0005807	10	171.5901	-322.7942	23.90753	6.383472e-06
4	-2.482182	-2.61e-05	-0.0016774	0.0018484	-1.6e-06	-0.0042020	0.0024819	NA	NA	1e-07	NA	NA	9	169.2317	-320.1481	26.55354	1.700135e-06
3	-2.481241	-2.64e-05	-0.0781279	0.0017918	-1.6e-06	-0.0040585	0.0028184	-1.98e-05	0.0030037	NA	NA	NA	10	167.3597	-314.3334	32.36824	9.286167e-08
9	-2.475019	-2.79e-05	0.0034649	0.0016366	-1.7e-06	-0.0041007	0.0025958	-8.70e-06	NA	NA	-0.0002856	NA	10	167.1645	-313.9430	32.75864	7.639465e-08
17	-2.474660	-2.78e-05	0.0034158	0.0016188	-1.7e-06	-0.0040574	0.0026121	-8.30e-06	NA	NA	NA	-0.0005125	10	166.6862	-312.9865	33.71517	4.735391e-08
26	-2.484329	-2.53e-05	-0.0752360	0.0019053	-1.6e-06	-0.0040689	0.0028130	NA	0.0026225	NA	-0.0021992	0.0000193	11	167.1639	-311.8639	34.83783	2.701296e-08
6	-2.491163	-2.51e-05	-0.0658585	0.0021010	-1.6e-06	-0.0040506	0.0025469	NA	0.0020406	7e-07	NA	NA	10	164.6871	-308.9882	37.71345	6.414165e-09
12	-2.477129	-2.70e-05	0.0008182	0.0017063	-1.7e-06	-0.0041453	0.0026049	NA	NA	-1e-07	-0.0012197	NA	10	164.3681	-308.3503	38.35143	4.662337e-09
20	-2.478881	-2.67e-05	-0.0012593	0.0017414	-1.6e-06	-0.0040546	0.0025662	NA	NA	1e-07	NA	-0.0009401	10	164.1511	-307.9162	38.78547	3.752787e-09
11	-2.481302	-2.63e-05	-0.0819535	0.0017988	-1.6e-06	-0.0040313	0.0028522	-1.59e-05	0.0030845	NA	-0.0007358	NA	11	162.8293	-303.1947	43.50702	3.540640e-10
19	-2.481434	-2.64e-05	-0.0793184	0.0017982	-1.6e-06	-0.0040698	0.0028161	-2.04e-05	0.0030529	NA	NA	0.0001175	11	162.3563	-302.2485	44.45314	2.206145e-10
25	-2.474659	-2.78e-05	0.0033284	0.0016192	-1.7e-06	-0.0040567	0.0026133	-8.10e-06	NA	NA	-0.0000445	-0.0005046	11	162.1775	-301.8910	44.81064	1.845031e-10
5	-2.474008	-2.75e-05	0.0074117	0.0016021	-1.7e-06	-0.0041954	0.0026138	-1.16e-05	NA	-3e-07	NA	NA	10	159.5612	-298.7363	47.96534	3.810396e-11
14	-2.486942	-2.61e-05	-0.0945190	0.0019825	-1.6e-06	-0.0038735	0.0027955	NA	0.0030906	6e-07	-0.0021341	NA	11	160.0554	-297.6469	49.05480	2.210023e-11
22	-2.487860	-2.57e-05	-0.0653996	0.0019940	-1.6e-06	-0.0039035	0.0026310	NA	0.0020393	7e-07	NA	-0.0009389	11	159.6068	-296.7496	49.95208	1.411094e-11
28	-2.476677	-2.70e-05	0.0003005	0.0016849	-1.7e-06	-0.0040728	0.0026183	NA	NA	0e+00	-0.0008419	-0.0005743	11	159.3941	-296.3242	50.37745	1.140740e-11
27	-2.481836	-2.63e-05	-0.0860582	0.0018177	-1.6e-06	-0.0040555	0.0028538	-1.66e-05	0.0032358	NA	-0.0009089	0.0003162	12	157.8956	-291.2420	55.45971	8.986463e-13
7	-2.483006	-2.68e-05	-0.0885776	0.0018448	-1.6e-06	-0.0039510	0.0028016	-1.91e-05	0.0032375	4e-07	NA	NA	11	155.2786	-288.0933	58.60841	1.861469e-13
13	-2.473997	-2.75e-05	0.0067813	0.0016039	-1.7e-06	-0.0041846	0.0026225	-1.01e-05	NA	-3e-07	-0.0002506	NA	11	155.0117	-287.5594	59.14227	1.425376e-13
21	-2.473936	-2.75e-05	0.0059245	0.0015963	-1.7e-06	-0.0041256	0.0026299	-9.40e-06	NA	-2e-07	NA	-0.0004528	11	154.5518	-286.6396	60.06208	8.999034e-14
30	-2.486908	-2.61e-05	-0.0943380	0.0019813	-1.6e-06	-0.0038720	0.0027955	NA	0.0030843	6e-07	-0.0021212	-0.0000167	12	155.0948	-285.6402	61.06147	5.459869e-14
15	-2.483310	-2.67e-05	-0.0942672	0.0018598	-1.6e-06	-0.0039062	0.0028373	-1.45e-05	0.0033594	4e-07	-0.0008255	NA	12	150.7592	-276.9690	69.73265	7.149220e-16
23	-2.483076	-2.68e-05	-0.0890223	0.0018472	-1.6e-06	-0.0039590	0.0028007	-1.94e-05	0.0032594	4e-07	NA	0.0000624	12	150.2835	-276.0176	70.68404	4.442896e-16
29	-2.473936	-2.75e-05	0.0058337	0.0015967	-1.7e-06	-0.0041249	0.0026312	-9.20e-06	NA	-2e-07	-0.0000470	-0.0004445	12	150.0439	-275.5386	71.16310	3.496553e-16
31	-2.483668	-2.67e-05	-0.0971660	0.0018729	-1.6e-06	-0.0039327	0.0028394	-1.52e-05	0.0034748	4e-07	-0.0009683	0.0002686	13	145.8290	-265.0161	81.68558	1.814318e-18

kable(Armenta.dredged.extinct.model, "html", caption = "BAMM Extinction Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM Extinction Spectrophotometry Dichromatism Dredge Table

	(Intercept)	bioclim4	Colour.discriminability	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	bioclim4:Colour.discriminability	Colour.discriminability:log(range.size.m2)	Colour.discriminability:NPP	Colour.discriminability:PC1.LIG	Colour.discriminability:residuals.PC1	df	logLik	AICc	delta	weight
0	-5.638235	-2.60e-06	-0.0007063	0.0098190	-7.0e-07	-0.0158136	0.0047533	NA	NA	NA	NA	NA	8	-295.8250	607.9017	0.000000	9.760815e-01
8	-5.603206	-9.50e-06	0.0013768	0.0086377	-1.0e-06	-0.0153132	0.0048284	NA	NA	NA	-0.0069542	NA	9	-299.3662	617.0476	9.145821	1.008083e-02
2	-5.660518	3.00e-06	-0.1231044	0.0104727	-5.0e-07	-0.0158420	0.0050480	NA	0.0042674	NA	NA	NA	9	-299.6673	617.6499	9.748116	7.459498e-03
16	-5.622155	-4.50e-06	0.0006290	0.0091959	-8.0e-07	-0.0155208	0.0048853	NA	NA	NA	NA	-0.0039732	9	-299.8737	618.0626	10.160901	6.068405e-03
10	-5.634016	-7.00e-07	-0.2647213	0.0094634	-8.0e-07	-0.0151204	0.0055100	NA	0.0093146	NA	-0.0104934	NA	10	-302.5940	625.5740	17.672233	1.419088e-04
1	-5.607926	-1.09e-05	0.0181308	0.0087106	-8.0e-07	-0.0155487	0.0048392	-4.15e-05	NA	NA	NA	NA	9	-304.4774	627.2700	19.368277	6.077406e-05
24	-5.601118	-9.30e-06	0.0016660	0.0085321	-1.0e-06	-0.0152558	0.0048753	NA	NA	NA	-0.0059789	-0.0017298	10	-303.5448	627.4756	19.573860	5.483734e-05
18	-5.643865	9.00e-07	-0.1162808	0.0098369	-6.0e-07	-0.0155557	0.0051632	NA	0.0040748	NA	NA	-0.0038665	10	-303.7310	627.8479	19.946169	4.552297e-05
4	-5.638385	1.60e-06	0.0142457	0.0097934	-7.0e-07	-0.0164385	0.0048215	NA	NA	-1.6e-06	NA	NA	9	-307.2012	632.7177	24.815946	3.988153e-06
3	-5.634989	-4.40e-06	-0.2351668	0.0093492	-6.0e-07	-0.0154226	0.0055422	-7.09e-05	0.0092968	NA	NA	NA	10	-307.7169	635.8199	27.918123	8.455567e-07
26	-5.634738	-6.00e-07	-0.2680032	0.0094910	-8.0e-07	-0.0151276	0.0055105	NA	0.0094278	NA	-0.0106989	0.0002882	11	-306.7717	636.0073	28.105536	7.699216e-07
9	-5.601830	-1.01e-05	0.0038916	0.0085790	-1.0e-06	-0.0153175	0.0048344	-5.90e-06	NA	NA	-0.0063678	NA	10	-307.9890	636.3640	28.462222	6.441591e-07
17	-5.605057	-1.01e-05	0.0144751	0.0085837	-9.0e-07	-0.0154279	0.0049016	-3.16e-05	NA	NA	NA	-0.0024957	10	-308.6217	637.6294	29.727669	3.421408e-07
12	-5.597342	-2.80e-06	0.0298211	0.0083824	-1.2e-06	-0.0163988	0.0049698	NA	NA	-3.1e-06	-0.0081742	NA	10	-310.5363	641.4586	33.556851	5.043227e-08
6	-5.658376	3.70e-06	-0.1065804	0.0104008	-5.0e-07	-0.0160276	0.0050397	NA	0.0038485	-5.0e-07	NA	NA	10	-310.9671	642.3201	34.418384	3.278151e-08
20	-5.621991	0.00e+00	0.0167326	0.0091558	-8.0e-07	-0.0161867	0.0049613	NA	NA	-1.8e-06	NA	-0.0040535	10	-311.2374	642.8608	34.959076	2.501608e-08
11	-5.629948	-2.80e-06	-0.2763034	0.0092454	-7.0e-07	-0.0151246	0.0056014	-2.96e-05	0.0101592	NA	-0.0078889	NA	11	-311.1065	644.6769	36.775141	1.008941e-08
19	-5.633651	-4.40e-06	-0.2294802	0.0093043	-6.0e-07	-0.0153988	0.0055382	-6.79e-05	0.0090580	NA	NA	-0.0005592	11	-311.9033	646.2705	38.368752	4.547968e-09
25	-5.600503	-9.60e-06	0.0028601	0.0085062	-1.0e-06	-0.0152591	0.0048772	-2.80e-06	NA	NA	-0.0057187	-0.0016952	11	-312.1582	646.7804	38.878672	3.524436e-09
14	-5.629632	7.00e-07	-0.2337836	0.0093153	-9.0e-07	-0.0154751	0.0054976	NA	0.0085422	-1.0e-06	-0.0105806	NA	11	-313.8788	650.2216	42.319887	6.307240e-10
5	-5.597683	-4.80e-06	0.0563682	0.0082702	-1.0e-06	-0.0167800	0.0050137	-5.60e-05	NA	-3.5e-06	NA	NA	10	-315.5658	651.5177	43.615917	3.299210e-10
28	-5.595913	-2.90e-06	0.0291288	0.0083088	-1.2e-06	-0.0163194	0.0050015	NA	NA	-3.0e-06	-0.0073794	-0.0013401	11	-314.7274	651.9188	44.017031	2.699661e-10
22	-5.640481	1.90e-06	-0.0917470	0.0097224	-6.0e-07	-0.0158266	0.0051523	NA	0.0034525	-7.0e-07	NA	-0.0039161	11	-315.0223	652.5086	44.606904	2.010110e-10
27	-5.631914	-2.80e-06	-0.2880311	0.0093148	-7.0e-07	-0.0151484	0.0056114	-3.23e-05	0.0106050	NA	-0.0083138	0.0009354	12	-315.2432	655.0356	47.133889	5.681881e-11
7	-5.626706	-2.30e-06	-0.1850096	0.0090594	-7.0e-07	-0.0160232	0.0055363	-7.41e-05	0.0081223	-1.6e-06	NA	NA	11	-318.9526	660.3691	52.467379	3.947659e-12
13	-5.591850	-4.10e-06	0.0418431	0.0081474	-1.2e-06	-0.0165353	0.0050064	-2.08e-05	NA	-3.4e-06	-0.0062482	NA	11	-319.0855	660.6349	52.733197	3.456350e-12
30	-5.630279	8.00e-07	-0.2368065	0.0093399	-9.0e-07	-0.0154785	0.0054981	NA	0.0086443	-1.0e-06	-0.0107548	0.0002452	12	-318.0553	660.6599	52.758122	3.413543e-12
21	-5.596145	-4.60e-06	0.0511733	0.0082013	-1.0e-06	-0.0166104	0.0050502	-4.75e-05	NA	-3.3e-06	NA	-0.0018999	11	-319.7398	661.9436	54.041824	1.796607e-12
15	-5.623288	-1.20e-06	-0.2342121	0.0090106	-8.0e-07	-0.0156238	0.0055951	-3.32e-05	0.0091760	-1.4e-06	-0.0076945	NA	12	-322.3572	669.2638	61.362018	4.622673e-14
23	-5.625559	-2.30e-06	-0.1802420	0.0090210	-7.0e-07	-0.0159975	0.0055326	-7.14e-05	0.0079155	-1.6e-06	NA	-0.0005036	12	-323.1385	670.8264	62.924614	2.116312e-14
29	-5.591361	-4.00e-06	0.0398539	0.0081165	-1.1e-06	-0.0164555	0.0050275	-1.84e-05	NA	-3.3e-06	-0.0058413	-0.0010734	12	-323.2685	671.0863	63.184602	1.858335e-14
31	-5.625254	-1.10e-06	-0.2459261	0.0090800	-8.0e-07	-0.0156498	0.0056052	-3.60e-05	0.0096237	-1.4e-06	-0.0081240	0.0009472	13	-326.4928	679.6277	71.725922	2.596570e-16

Run the top model and 100 models for each one:

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
Armenta.MCC.DR.top <- gls(MCC.DR ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.DR.top, 'data/Armenta.MCC.DR.top.rds')

#Run model for ND
Armenta.MCC.ND.top <- gls(MCC.ND ~ Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.ND.top, 'data/Armenta.MCC.ND.top.rds')

#Run the 100 models for DR and ND using the best model:
Armenta.data.noMCC <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), Armenta.data, by = "binomial") %>% filter(Colour.discriminability != "NA")
#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Armenta.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.DR.model.data <- lapply(Armenta.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Armenta.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(MCC.Armenta.model.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Armenta.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(1, phy = y, fixed = FALSE), 
    data = x, 
    method = "REML")
}, x = Armenta.DR.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.DR.pgls.models, "data/Armenta.DR.pgls.models.rds")

#Now for Node Density:
Armenta.ND.model.data <- lapply(nd.list, function(x) {
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.ND.model.data <- lapply(Armenta.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Armenta.ND.model.data, y = Armenta.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.ND.pgls.models, "data/Armenta.ND.pgls.models.rds")

#Run the BAMM models

Armenta.MCC.Lambda.top <- gls(mean.lambda ~ Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Lambda.top, 'data/Armenta.MCC.Lambda.top.rds')

Armenta.MCC.Mu.top <- gls(mean.mu ~ Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Mu.top, 'data/Armenta.MCC.Mu.top.rds')

Armenta.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.BAMM.model.data <- lapply(Armenta.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.lambda.pgls.models, "data/Armenta.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Armenta.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ Colour.discriminability
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.mu.pgls.models, "data/Armenta.BAMM.mu.pgls.models.rds")

Using the MCC models as well as the 100 PGLS models we can generate the plots similar to Supp Fig S8:

Armenta.DR.pgls.models <- readRDS("data/Armenta.DR.pgls.models.rds")
Armenta.MCC.DR.top <- readRDS('data/Armenta.MCC.DR.top.rds')

Armenta.DR.pgls.summary <- lapply(Armenta.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.DR.summary <- data.frame(Armenta.MCC.DR.top$coefficients, confint(Armenta.MCC.DR.top)) %>% tibble::rownames_to_column()

Armenta.DR.pgls.summary <- bind_rows(Armenta.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `Colour.discriminability` = "Sexual Dichromatism"
                    )

Armenta.DR.pgls.summary$Parameter = factor(Armenta.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.DR.summary$Parameter = factor(Armenta.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.DR.plot <-Armenta.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.DR.plot <- Armenta.DR.plot + geom_errorbarh(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta DR")

####################################
#For ND

Armenta.ND.pgls.models <- readRDS("data/Armenta.ND.pgls.models.rds")
Armenta.MCC.ND.top <- readRDS('data/Armenta.MCC.ND.top.rds')

Armenta.ND.pgls.summary <- lapply(Armenta.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.ND.summary <- data.frame(Armenta.MCC.ND.top$coefficients, confint(Armenta.MCC.ND.top)) %>% tibble::rownames_to_column()

Armenta.ND.pgls.summary <- bind_rows(Armenta.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.ND.pgls.summary$Parameter = factor(Armenta.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.ND.summary$Parameter = factor(Armenta.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.ND.plot <-Armenta.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.ND.plot <- Armenta.ND.plot + geom_errorbarh(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta ND")

######################

#For Lambda
Armenta.BAMM.lambda.pgls.models <- readRDS("data/Armenta.BAMM.lambda.pgls.models.rds")
Armenta.MCC.Lambda.top <- readRDS('data/Armenta.MCC.Lambda.top.rds')

Armenta.Lambda.pgls.summary <- lapply(Armenta.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Lambda.summary <- data.frame(Armenta.MCC.Lambda.top$coefficients, confint(Armenta.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Armenta.Lambda.pgls.summary <- bind_rows(Armenta.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Lambda.pgls.summary$Parameter = factor(Armenta.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.Lambda.summary$Parameter = factor(Armenta.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.Lambda.plot <-Armenta.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Lambda.plot <- Armenta.Lambda.plot + geom_errorbarh(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Speciaton")

######################

#For Mu

Armenta.BAMM.Mu.pgls.models <- readRDS("data/Armenta.BAMM.mu.pgls.models.rds")
Armenta.MCC.Mu.top <- readRDS('data/Armenta.MCC.Mu.top.rds')

Armenta.Mu.pgls.summary <- lapply(Armenta.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Mu.summary <- data.frame(Armenta.MCC.Mu.top$coefficients, confint(Armenta.MCC.Mu.top)) %>% tibble::rownames_to_column()

Armenta.Mu.pgls.summary <- bind_rows(Armenta.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Mu.pgls.summary$Parameter = factor(Armenta.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.MCC.Mu.summary$Parameter = factor(Armenta.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Colour.discriminability'))

Armenta.Mu.plot <-Armenta.Mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Mu.plot <- Armenta.Mu.plot + geom_errorbarh(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Extinction")

Plot the results

grid.arrange(symmetrise_scale(Armenta.DR.plot, "x"),
             symmetrise_scale(Armenta.ND.plot, "x"),
             symmetrise_scale(Armenta.Lambda.plot, "x"), 
             symmetrise_scale(Armenta.Mu.plot, "x"), 
             nrow = 1)

Figure S10: Measures of sexual dichromatism using a more precise measure in spectrophotometry do not change the main patterns seen within the PGLS models analysed here. This dataset is a subset of the complete dataset (n = 581), thus drawing conclusions for the other predictors (e.g. borderline long term temperature variation and spatial temperature variation) potentially risks type I error; additionally these environmental predictors were not measured differently in this analysis compared to the previous above so we have no obvious reason to make conclusions from this subset for the effects of the environmental predictors.

Table S11: MCC estimates from the above plot are presented as numerical values below.

Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates from Spec measures")

MCC DR Estimates from Spec measures

Parameter	Estimate	LCI	UCI
Colour.discriminability	0.007071	-0.0482	0.06234
log(range.size.m2)	0.009034	-0.02176	0.03982
bioclim4	0.0001789	-8.595e-05	0.0004437
residuals.PC1	0.0156	-0.0073	0.0385
PC1.LIG	-0.03157	-0.07051	0.007379
NPP	7.814e-06	-5.345e-06	2.097e-05

Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates from Spec measures")

MCC ND Estimates from Spec measures

Parameter	Estimate	LCI	UCI
Colour.discriminability	8.285e-05	-0.003516	0.003682
log(range.size.m2)	-3.776e-05	-0.001739	0.001663
bioclim4	6.941e-06	-8.315e-06	2.22e-05
residuals.PC1	0.001053	-0.0001951	0.002302
PC1.LIG	-0.0008841	-0.003057	0.001288
NPP	4.626e-07	-2.874e-07	1.213e-06

Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates from Spec measures")

MCC BAMM Speciation Estimates from Spec measures

Parameter	Estimate	LCI	UCI
Colour.discriminability	-0.0005285	-0.01583	0.01478
log(range.size.m2)	0.001845	-0.005712	0.009402
bioclim4	-2.585e-05	-9.338e-05	4.167e-05
residuals.PC1	0.002489	-0.003265	0.008243
PC1.LIG	-0.004244	-0.01356	0.005072
NPP	-1.636e-06	-4.926e-06	1.654e-06

Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates from Spec measures")

MCC BAMM Extinction Estimates from Spec measures

Parameter	Estimate	LCI	UCI
Colour.discriminability	-0.0007063	-0.03742	0.036
log(range.size.m2)	0.009819	-0.007562	0.0272
bioclim4	-2.645e-06	-0.0001577	0.0001524
residuals.PC1	0.004753	-0.008041	0.01755
PC1.LIG	-0.01581	-0.03777	0.006146
NPP	-6.949e-07	-8.335e-06	6.946e-06

Table S12: Highest posterior density intervals are calculated for the above figure (ie; models using Spectrophotometry values of sexual dichromatism). The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -). These intervals are calculated in the same way as in Table S7.

Armenta.hpd.DR.top <- list()
Armenta.DR.pgls.summary <- na.omit(Armenta.DR.pgls.summary)
for(x in unique(Armenta.DR.pgls.summary$Parameter)) {
Armenta.hpd.DR.top[[x]] = hdi(Armenta.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.DR.top <- bind_rows(Armenta.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Armenta.hpd.DR.top, 'data/Armenta.hpd.DR.top.rds')

#For ND
Armenta.hpd.ND.top <- list()
Armenta.ND.pgls.summary <- na.omit(Armenta.ND.pgls.summary)
for(x in unique(Armenta.ND.pgls.summary$Parameter)) {
Armenta.hpd.ND.top[[x]] = hdi(Armenta.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.ND.top <- bind_rows(Armenta.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.ND.top, 'data/Armenta.hpd.ND.top.rds')

###############
Armenta.hpd.Lambda.top <- list()
Armenta.Lambda.pgls.summary <- na.omit(Armenta.Lambda.pgls.summary)
for(x in unique(Armenta.Lambda.pgls.summary$Parameter)) {
Armenta.hpd.Lambda.top[[x]] = hdi(Armenta.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Lambda.top <- bind_rows(Armenta.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Armenta.hpd.Lambda.top, 'data/Armenta.hpd.Lambda.top.rds')

################
Armenta.hpd.Mu.top <- list()
Armenta.Mu.pgls.summary <- na.omit(Armenta.Mu.pgls.summary)
for(x in unique(Armenta.Mu.pgls.summary$Parameter)) {
Armenta.hpd.Mu.top[[x]] = hdi(Armenta.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Mu.top <- bind_rows(Armenta.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.Mu.top, 'data/Armenta.hpd.Mu.top.rds')

Armenta.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta DR HPD Interval")

Armenta DR HPD Interval

	Colour.discriminability	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0178	-0.0148	-0.000139	0.0108	-0.0424	-4.23e-06
Upper	0.0349	0.0316	0.00016	0.0347	0.00792	1.45e-05

Armenta.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta ND HPD Interval")

Armenta ND HPD Interval

	Colour.discriminability	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.00176	-0.000596	-4.97e-06	5e-04	-0.0023	-7e-08
Upper	0.00173	0.00133	1.37e-05	0.00184	0.000556	6.84e-07

Armenta.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Speciation HPD Interval")

Armenta BAMM Speciation HPD Interval

	Colour.discriminability	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0127	-0.00713	-4.41e-05	-0.00289	-0.0173	-4.05e-06
Upper	0.011	0.00837	5.98e-05	0.00677	0.00505	3.69e-06

Armenta.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Extinction HPD Interval")

Armenta BAMM Extinction HPD Interval

	Colour.discriminability	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0282	-0.019	-0.000119	-0.00577	-0.0377	-1.05e-05
Upper	0.0172	0.0207	0.000136	0.0158	0.00783	1.37e-05

Analysis using male-bias measure of sexual selection

Sexual dichromatism is expected to be a good measure of sexual selection strength in birds. However the relationship between sexual dichromatism and male-biased measures of sexual selection (social mating system, sexual size dimorphism, and paternal care) is expected to be relatively noisy given the precision of the measurment used (Dale et al. 2015). Here we use a dataset of sexual selection for 2,465 species within the orginal dataset of sexual dichromatism scores from RGB measures. This male-bias sexual selection score is based on three components, combined in a phylogenetic pca (ppca) with the following loadings:

• sexual size dimorphism: 0.37
• social polygyny: 0.57
• paternal care: -0.57

As such, high values for this score are expected to have high sexual selection (high dimorphism, high polygyny and low paternal care). PGLS models using this dataset are calculated in the same process as above.

SS.subset <- plumage.scores %>% filter(Sexual_selection_ppca != "NA")
SS.subset <- inner_join(SS.subset, restricted.data %>% dplyr::select(TipLabel, binomial, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), by = "TipLabel")
SS.subset <- inner_join(SS.subset, MCC.BAMM.df, by = "TipLabel")

SS.subset %>% ggplot(aes(x = SDi, y = Sexual_selection_ppca))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_minimal()

Figure S11: The relationship between absolute sexual dichromatism and the sexual selection. Although positively correlated the distribution is very noisy with species having high sexual dichromatism and low sexual selection and visa versa.Colour.discriminability

# summary(lm(SDi ~ Sexual_selection_ppcaColour.discriminability
#    data = SS.subset))$r.squared %>% pander()

cor(SS.subset$SDi, SS.subset$Sexual_selection_ppca)

## [1] 0.3414535

#Prune tree
pruned.MCC.Subset.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(SS.subset$TipLabel, MCC.passerine$tip.label)])

saveRDS(pruned.MCC.Subset.tree, 'data/pruned.MCC.Subset.tree.rds')

#Set rownames to match tree
rownames(SS.subset) <- SS.subset$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Subset.global <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Subset.global <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Subset.global <- gls(mean.lambda ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variatioColour.discriminability
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Subset.global <- gls(mean.mu ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ sqrt(var.mu),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("SS.subset"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Subset.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Subset.dredged.ND.model <- pdredge(MCC.ND.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.ND.model, 'data/Subset.dredged.ND.model.rds')

Subset.dredged.DR.model <- pdredge(MCC.DR.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.DR.model, 'data/Subset.dredged.DR.model.rds')

Subset.dredged.spec.model <- pdredge(MCC.Lambda.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.spec.model.rds')

Subset.dredged.extinct.model <- pdredge(MCC.Extinction.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.extinct.model.rds')

Table S13: All top models (\(\delta AICc > 4\)) using male-biased sexual selection measures do not contain interactions between sexual selection and the environmental variables. Thus interaction terms are not included in further analysis.

Subset.dredged.ND.model <- readRDS("data/Subset.dredged.ND.model.rds")
Subset.dredged.DR.model <- readRDS("data/Subset.dredged.DR.model.rds")
Subset.dredged.spec.model <- readRDS("data/Subset.dredged.spec.model.rds")
Subset.dredged.extinct.model <- readRDS("data/Subset.dredged.extinct.model.rds")


kable(Subset.dredged.ND.model, "html", caption = "ND Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

ND Male-bias Sexual Selection Dredge Table

	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	Sexual_selection_ppca	bioclim4:Sexual_selection_ppca	log(range.size.m2):Sexual_selection_ppca	NPP:Sexual_selection_ppca	PC1.LIG:Sexual_selection_ppca	residuals.PC1:Sexual_selection_ppca	df	logLik	AICc	delta	weight
0	0.1484807	-5e-07	-0.0001800	0	0.0001073	0.0000788	0.0004381	NA	NA	NA	NA	NA	8	5451.548	-10887.04	0.00000	9.994587e-01
16	0.1482246	-5e-07	-0.0001720	0	0.0000929	0.0000644	0.0003782	NA	NA	NA	NA	-0.0001167	9	5444.064	-10870.05	16.98341	2.050523e-04
8	0.1479484	-5e-07	-0.0001570	0	0.0000777	0.0001046	0.0004928	NA	NA	NA	0.0001280	NA	9	5444.026	-10869.98	17.05889	1.974574e-04
2	0.1483243	-5e-07	-0.0001765	0	0.0001017	0.0000834	0.0018018	NA	-4.87e-05	NA	NA	NA	9	5443.645	-10869.22	17.82105	1.348878e-04
1	0.1478764	-9e-07	-0.0001561	0	0.0000738	0.0000871	0.0012068	-2.3e-06	NA	NA	NA	NA	9	5440.065	-10862.06	24.98016	3.761776e-06
4	0.1484220	-6e-07	-0.0001730	0	0.0001024	0.0000814	0.0008188	NA	NA	0e+00	NA	NA	9	5436.296	-10854.52	32.51872	8.677880e-08
24	0.1477981	-5e-07	-0.0001528	0	0.0000692	0.0000925	0.0004458	NA	NA	NA	0.0001196	-0.0000845	10	5436.455	-10852.82	34.21679	3.712640e-08
18	0.1479767	-6e-07	-0.0001662	0	0.0000832	0.0000688	0.0022340	NA	-6.65e-05	NA	NA	-0.0001323	10	5436.226	-10852.36	34.67399	2.953948e-08
10	0.1478037	-5e-07	-0.0001538	0	0.0000725	0.0001088	0.0017783	NA	-4.59e-05	NA	0.0001274	NA	10	5436.116	-10852.14	34.89512	2.644759e-08
17	0.1478092	-9e-07	-0.0001543	0	0.0000700	0.0000806	0.0011337	-2.2e-06	NA	NA	NA	-0.0000480	10	5432.459	-10844.83	42.20969	6.824126e-10
9	0.1474813	-8e-07	-0.0001389	0	0.0000518	0.0001085	0.0011728	-2.1e-06	NA	NA	0.0001104	NA	10	5432.373	-10844.66	42.38150	6.262388e-10
3	0.1478909	-9e-07	-0.0001563	0	0.0000742	0.0000866	0.0010291	-2.3e-06	6.70e-06	NA	NA	NA	10	5432.150	-10844.21	42.82661	5.012840e-10
20	0.1482076	-6e-07	-0.0001673	0	0.0000906	0.0000682	0.0006967	NA	NA	0e+00	NA	-0.0001026	10	5428.776	-10837.46	49.57473	1.716915e-11
12	0.1478969	-6e-07	-0.0001505	0	0.0000733	0.0001068	0.0008569	NA	NA	0e+00	0.0001269	NA	10	5428.762	-10837.43	49.60299	1.692825e-11
6	0.1482266	-6e-07	-0.0001680	0	0.0000951	0.0000872	0.0025093	NA	-5.90e-05	0e+00	NA	NA	10	5428.428	-10836.77	50.27032	1.212563e-11
26	0.1475846	-6e-07	-0.0001479	0	0.0000610	0.0000958	0.0021036	NA	-5.95e-05	NA	0.0001173	-0.0000990	11	5428.596	-10835.09	51.95195	5.230492e-12
5	0.1477933	-1e-06	-0.0001475	0	0.0000674	0.0000903	0.0016611	-2.4e-06	NA	-1e-07	NA	NA	10	5424.857	-10829.62	57.41215	3.410950e-13
25	0.1474541	-8e-07	-0.0001383	0	0.0000503	0.0001049	0.0011368	-2.0e-06	NA	NA	0.0001086	-0.0000240	11	5424.751	-10827.39	59.64290	1.118085e-13
19	0.1477962	-9e-07	-0.0001541	0	0.0000696	0.0000808	0.0012568	-2.1e-06	-4.70e-06	NA	NA	-0.0000499	11	5424.591	-10827.07	59.96353	9.524668e-14
11	0.1474880	-8e-07	-0.0001390	0	0.0000520	0.0001082	0.0010942	-2.1e-06	2.90e-06	NA	0.0001103	NA	11	5424.457	-10826.81	60.23059	8.334079e-14
28	0.1477791	-6e-07	-0.0001479	0	0.0000668	0.0000965	0.0007717	NA	NA	0e+00	0.0001201	-0.0000699	11	5421.172	-10820.24	66.80118	3.119414e-15
22	0.1479328	-7e-07	-0.0001603	0	0.0000797	0.0000735	0.0027753	NA	-7.32e-05	0e+00	NA	-0.0001181	11	5420.966	-10819.82	67.21289	2.539046e-15
14	0.1477158	-6e-07	-0.0001460	0	0.0000665	0.0001122	0.0024557	NA	-5.58e-05	0e+00	0.0001260	NA	11	5420.885	-10819.66	67.37476	2.341650e-15
21	0.1477617	-1e-06	-0.0001470	0	0.0000656	0.0000868	0.0016033	-2.3e-06	NA	0e+00	NA	-0.0000251	11	5417.249	-10812.39	74.64696	6.171400e-17
13	0.1474084	-1e-06	-0.0001309	0	0.0000460	0.0001112	0.0016076	-2.1e-06	NA	0e+00	0.0001085	NA	11	5417.149	-10812.19	74.84702	5.583943e-17
7	0.1477851	-1e-06	-0.0001473	0	0.0000671	0.0000906	0.0017584	-2.4e-06	-3.60e-06	-1e-07	NA	NA	11	5416.948	-10811.79	75.24841	4.568560e-17
27	0.1474468	-8e-07	-0.0001381	0	0.0000501	0.0001050	0.0012077	-2.0e-06	-2.70e-06	NA	0.0001085	-0.0000251	12	5416.883	-10809.64	77.39880	1.558939e-17
30	0.1475395	-6e-07	-0.0001419	0	0.0000575	0.0001007	0.0026498	NA	-6.62e-05	0e+00	0.0001176	-0.0000846	12	5413.338	-10802.55	84.48797	4.502300e-19
29	0.1474072	-1e-06	-0.0001309	0	0.0000459	0.0001110	0.0016049	-2.1e-06	NA	0e+00	0.0001084	-0.0000012	12	5409.540	-10794.95	92.08377	1.009318e-20
23	0.1477348	-1e-06	-0.0001465	0	0.0000647	0.0000871	0.0018596	-2.2e-06	-9.70e-06	0e+00	NA	-0.0000288	12	5409.384	-10794.64	92.39589	8.634816e-21
15	0.1473924	-1e-06	-0.0001306	0	0.0000455	0.0001118	0.0017935	-2.1e-06	-6.80e-06	0e+00	0.0001087	NA	12	5409.241	-10794.35	92.68281	7.480788e-21
31	0.1473862	-1e-06	-0.0001305	0	0.0000452	0.0001112	0.0018080	-2.1e-06	-7.70e-06	0e+00	0.0001084	-0.0000042	13	5401.675	-10777.20	109.83562	1.410139e-24

kable(Subset.dredged.DR.model, "html", caption = "DR Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

DR Male-bias Sexual Selection Dredge Table

	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	Sexual_selection_ppca	bioclim4:Sexual_selection_ppca	log(range.size.m2):Sexual_selection_ppca	NPP:Sexual_selection_ppca	PC1.LIG:Sexual_selection_ppca	residuals.PC1:Sexual_selection_ppca	df	logLik	AICc	delta	weight
0	-2.644369	2.45e-05	-0.0113055	-7.0e-07	0.0068813	0.0004409	0.0388722	NA	NA	NA	NA	NA	9	-1717.437	3452.947	0.00000	9.851357e-01
16	-2.639661	2.53e-05	-0.0115625	-7.0e-07	0.0066554	-0.0000035	0.0381221	NA	NA	NA	NA	-0.0043664	10	-1721.450	3462.989	10.04197	6.499963e-03
2	-2.651446	2.27e-05	-0.0111253	-7.0e-07	0.0068731	0.0003113	0.1032726	NA	-0.0023191	NA	NA	NA	10	-1721.901	3463.892	10.94493	4.138415e-03
8	-2.643685	2.47e-05	-0.0113672	-8.0e-07	0.0067933	0.0002826	0.0391218	NA	NA	NA	-0.0017643	NA	10	-1721.915	3463.919	10.97215	4.082471e-03
1	-2.644777	1.73e-05	-0.0112474	-8.0e-07	0.0069016	0.0002884	0.0537640	-4.40e-05	NA	NA	NA	NA	10	-1726.048	3472.186	19.23901	6.543319e-05
18	-2.648357	2.30e-05	-0.0113613	-7.0e-07	0.0066118	-0.0002393	0.1227636	NA	-0.0030524	NA	NA	-0.0050116	11	-1725.773	3473.654	20.70712	3.140527e-05
24	-2.639316	2.54e-05	-0.0115972	-7.0e-07	0.0066022	-0.0000984	0.0383393	NA	NA	NA	-0.0012567	-0.0041819	11	-1725.963	3474.034	21.08668	2.597662e-05
10	-2.650479	2.30e-05	-0.0111903	-7.0e-07	0.0067936	0.0001735	0.1006675	NA	-0.0022170	NA	-0.0015983	NA	11	-1726.392	3474.892	21.94541	1.690876e-05
4	-2.639136	2.37e-05	-0.0113014	-1.1e-06	0.0067364	0.0005491	0.0550570	NA	NA	-1.6e-06	NA	NA	10	-1729.238	3478.565	25.61792	2.695472e-06
17	-2.641237	1.95e-05	-0.0114476	-8.0e-07	0.0067343	0.0000040	0.0497120	-3.36e-05	NA	NA	NA	-0.0031418	11	-1730.243	3482.593	29.64645	3.596263e-07
9	-2.644302	1.77e-05	-0.0112841	-8.0e-07	0.0068566	0.0002153	0.0532073	-4.20e-05	NA	NA	-0.0008810	NA	11	-1730.569	3483.246	30.29926	2.594752e-07
3	-2.648526	1.70e-05	-0.0111570	-7.0e-07	0.0068950	0.0002315	0.0876844	-3.97e-05	-0.0012746	NA	NA	NA	11	-1730.598	3483.303	30.35634	2.521738e-07
26	-2.647793	2.32e-05	-0.0113950	-7.0e-07	0.0065727	-0.0003052	0.1206409	NA	-0.0029697	NA	-0.0009490	-0.0048550	12	-1730.301	3484.728	31.78144	1.236643e-07
20	-2.635570	2.45e-05	-0.0115406	-1.0e-06	0.0065504	0.0001210	0.0518306	NA	NA	-1.4e-06	NA	-0.0040359	11	-1733.315	3488.737	35.78960	1.666795e-08
6	-2.646854	2.14e-05	-0.0110830	-1.1e-06	0.0067038	0.0004113	0.1349625	NA	-0.0027849	-1.9e-06	NA	NA	11	-1733.615	3489.338	36.39048	1.234251e-08
12	-2.638336	2.39e-05	-0.0113656	-1.1e-06	0.0066414	0.0003852	0.0556198	NA	NA	-1.7e-06	-0.0018515	NA	11	-1733.706	3489.519	36.57214	1.127085e-08
19	-2.647224	1.97e-05	-0.0113350	-7.0e-07	0.0066756	-0.0001778	0.1100270	-2.27e-05	-0.0023104	NA	NA	-0.0040281	12	-1734.642	3493.410	40.46326	1.610686e-09
25	-2.640864	1.99e-05	-0.0114754	-8.0e-07	0.0067004	-0.0000514	0.0493311	-3.21e-05	NA	NA	-0.0007139	-0.0030923	12	-1734.769	3493.665	40.71787	1.418151e-09
11	-2.648042	1.75e-05	-0.0111939	-8.0e-07	0.0068500	0.0001586	0.0870861	-3.76e-05	-0.0012730	NA	-0.0008797	NA	12	-1735.118	3494.364	41.41693	9.998225e-10
5	-2.639640	1.65e-05	-0.0112442	-1.1e-06	0.0067596	0.0003957	0.0694708	-4.36e-05	NA	-1.6e-06	NA	NA	11	-1737.859	3497.825	44.87789	1.771687e-10
22	-2.644516	2.18e-05	-0.0113100	-1.1e-06	0.0064809	-0.0001171	0.1492103	NA	-0.0034126	-1.7e-06	NA	-0.0046899	12	-1737.559	3499.245	46.29775	8.710986e-11
28	-2.635111	2.47e-05	-0.0115772	-1.1e-06	0.0064896	0.0000208	0.0524230	NA	NA	-1.4e-06	-0.0013747	-0.0038249	12	-1737.818	3499.763	46.81566	6.723639e-11
14	-2.645817	2.17e-05	-0.0111498	-1.2e-06	0.0066197	0.0002691	0.1325514	NA	-0.0026832	-1.9e-06	-0.0016616	NA	12	-1738.100	3500.327	47.37978	5.071161e-11
27	-2.646817	2.00e-05	-0.0113622	-7.0e-07	0.0066447	-0.0002281	0.1093216	-2.14e-05	-0.0022967	NA	-0.0006587	-0.0039771	13	-1739.168	3504.485	51.53820	6.340417e-12
21	-2.637030	1.86e-05	-0.0114218	-1.1e-06	0.0066273	0.0001336	0.0642539	-3.45e-05	NA	-1.4e-06	NA	-0.0027665	12	-1742.094	3508.315	55.36799	9.343053e-13
7	-2.644376	1.61e-05	-0.0111162	-1.1e-06	0.0067357	0.0003290	0.1181574	-3.75e-05	-0.0017672	-1.8e-06	NA	NA	12	-1742.351	3508.828	55.88132	7.228072e-13
13	-2.639083	1.70e-05	-0.0112841	-1.1e-06	0.0067082	0.0003160	0.0690252	-4.13e-05	NA	-1.6e-06	-0.0009791	NA	12	-1742.374	3508.875	55.92791	7.061633e-13
30	-2.643870	2.20e-05	-0.0113459	-1.1e-06	0.0064361	-0.0001880	0.1471758	NA	-0.0033252	-1.7e-06	-0.0010505	-0.0045120	13	-1742.080	3510.308	57.36050	3.450013e-13
23	-2.643470	1.87e-05	-0.0112857	-1.1e-06	0.0065443	-0.0000596	0.1365454	-2.19e-05	-0.0026920	-1.6e-06	NA	-0.0037463	13	-1746.434	3519.016	66.06877	4.434500e-15
29	-2.636582	1.90e-05	-0.0114524	-1.1e-06	0.0065870	0.0000719	0.0640021	-3.28e-05	NA	-1.4e-06	-0.0008236	-0.0027038	13	-1746.614	3519.376	66.42937	3.702894e-15
15	-2.643820	1.66e-05	-0.0111564	-1.1e-06	0.0066838	0.0002487	0.1177972	-3.52e-05	-0.0017704	-1.8e-06	-0.0009854	NA	13	-1746.865	3519.878	66.93090	2.881620e-15
31	-2.642995	1.90e-05	-0.0113159	-1.1e-06	0.0065069	-0.0001168	0.1359737	-2.04e-05	-0.0026794	-1.6e-06	-0.0007727	-0.0036828	14	-1750.955	3530.081	77.13422	1.753934e-17

kable(Subset.dredged.spec.model, "html", caption = "BAMM Speciation Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM Speciation Male-bias Sexual Selection Dredge Table

	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	Sexual_selection_ppca	bioclim4:Sexual_selection_ppca	log(range.size.m2):Sexual_selection_ppca	NPP:Sexual_selection_ppca	PC1.LIG:Sexual_selection_ppca	residuals.PC1:Sexual_selection_ppca	df	logLik	AICc	delta	weight
0	-2.291100	-7.6e-06	0.0005471	-1e-07	-0.0005670	0.0004011	0.0009422	NA	NA	NA	NA	NA	8	2239.623	-4463.188	0.00000	9.982421e-01
16	-2.291695	-7.7e-06	0.0005655	-1e-07	-0.0005911	0.0003593	0.0008204	NA	NA	NA	NA	-0.0002658	9	2233.363	-4448.652	14.53554	6.964394e-04
8	-2.292345	-7.7e-06	0.0005965	-1e-07	-0.0006020	0.0004437	0.0009637	NA	NA	NA	0.0003152	NA	9	2233.113	-4448.152	15.03575	5.423296e-04
2	-2.290976	-7.6e-06	0.0005418	-1e-07	-0.0005524	0.0003885	-0.0024843	NA	0.0001231	NA	NA	NA	9	2233.054	-4448.035	15.15316	5.114101e-04
1	-2.291399	-7.9e-06	0.0005605	-1e-07	-0.0005809	0.0004077	0.0013441	-1.2e-06	NA	NA	NA	NA	9	2228.666	-4439.258	23.92997	6.351986e-06
24	-2.292817	-7.8e-06	0.0006107	-1e-07	-0.0006217	0.0004057	0.0008568	NA	NA	NA	0.0003034	-0.0002315	10	2226.838	-4433.586	29.60133	3.727258e-07
18	-2.291559	-7.6e-06	0.0005602	-1e-07	-0.0005786	0.0003533	-0.0016797	NA	0.0000902	NA	NA	-0.0002452	10	2226.797	-4433.504	29.68333	3.577527e-07
4	-2.291653	-7.9e-06	0.0005740	-1e-07	-0.0005907	0.0004065	0.0023362	NA	NA	-2e-07	NA	NA	9	2225.687	-4433.301	29.88723	3.230772e-07
10	-2.292220	-7.7e-06	0.0005912	-1e-07	-0.0005877	0.0004314	-0.0023393	NA	0.0001187	NA	0.0003140	NA	10	2226.542	-4432.994	30.19380	2.771635e-07
17	-2.291756	-7.8e-06	0.0005687	-1e-07	-0.0005941	0.0003635	0.0009507	-4.0e-07	NA	NA	NA	-0.0002518	10	2222.457	-4424.824	38.36346	4.663555e-09
3	-2.291393	-7.9e-06	0.0005605	-1e-07	-0.0005682	0.0003937	-0.0032200	-1.9e-06	0.0001720	NA	NA	NA	10	2222.169	-4424.249	38.93904	3.497280e-09
9	-2.292576	-7.9e-06	0.0006068	-1e-07	-0.0006129	0.0004487	0.0012922	-1.0e-06	NA	NA	0.0003117	NA	10	2222.151	-4424.212	38.97548	3.434136e-09
20	-2.292081	-7.9e-06	0.0005863	-1e-07	-0.0006079	0.0003718	0.0020931	NA	NA	-1e-07	NA	-0.0002167	10	2219.415	-4418.741	44.44700	2.226858e-10
26	-2.292680	-7.7e-06	0.0006054	-1e-07	-0.0006090	0.0003997	-0.0016523	NA	0.0000905	NA	0.0003035	-0.0002108	11	2220.273	-4418.437	44.75033	1.913486e-10
12	-2.292870	-8.0e-06	0.0006222	-1e-07	-0.0006246	0.0004485	0.0023252	NA	NA	-1e-07	0.0003115	NA	10	2219.171	-4418.253	44.93462	1.745048e-10
6	-2.291537	-7.8e-06	0.0005690	-1e-07	-0.0005789	0.0003970	-0.0002205	NA	0.0000895	-1e-07	NA	NA	10	2219.115	-4418.140	45.04770	1.649120e-10
19	-2.291677	-7.8e-06	0.0005670	-1e-07	-0.0005823	0.0003628	-0.0022333	-1.1e-06	0.0001231	NA	NA	-0.0001991	11	2215.992	-4409.877	53.31107	2.647794e-12
25	-2.292858	-7.8e-06	0.0006128	-1e-07	-0.0006237	0.0004086	0.0009470	-3.0e-07	NA	NA	0.0003030	-0.0002218	11	2215.932	-4409.756	53.43142	2.493164e-12
5	-2.292014	-8.2e-06	0.0005902	-1e-07	-0.0006074	0.0004143	0.0028381	-1.4e-06	NA	-2e-07	NA	NA	10	2214.737	-4409.384	53.80399	2.069413e-12
11	-2.292555	-8.0e-06	0.0006063	-1e-07	-0.0006006	0.0004350	-0.0029759	-1.7e-06	0.0001609	NA	0.0003077	NA	11	2215.650	-4409.191	53.99633	1.879667e-12
28	-2.293197	-8.0e-06	0.0006311	-1e-07	-0.0006382	0.0004180	0.0021203	NA	NA	-1e-07	0.0003025	-0.0001829	11	2212.889	-4403.671	59.51672	1.189444e-13
22	-2.291971	-7.9e-06	0.0005818	-1e-07	-0.0005983	0.0003671	0.0002469	NA	0.0000652	-1e-07	NA	-0.0002034	11	2212.849	-4403.591	59.59667	1.142835e-13
14	-2.292756	-7.9e-06	0.0006172	-1e-07	-0.0006133	0.0004393	-0.0001273	NA	0.0000859	-1e-07	0.0003107	NA	11	2212.598	-4403.089	60.09851	8.892208e-14
21	-2.292218	-8.1e-06	0.0005934	-1e-07	-0.0006146	0.0003809	0.0024002	-8.0e-07	NA	-1e-07	NA	-0.0001865	11	2208.522	-4394.936	68.25187	1.508452e-15
27	-2.292778	-7.9e-06	0.0006110	-1e-07	-0.0006121	0.0004077	-0.0021330	-9.0e-07	0.0001191	NA	0.0003020	-0.0001710	12	2209.466	-4394.805	68.38296	1.412746e-15
7	-2.291971	-8.2e-06	0.0005884	-1e-07	-0.0005954	0.0004025	-0.0009542	-2.0e-06	0.0001395	-1e-07	NA	NA	11	2208.232	-4394.357	68.83056	1.129455e-15
13	-2.293157	-8.2e-06	0.0006351	-1e-07	-0.0006382	0.0004545	0.0027477	-1.2e-06	NA	-2e-07	0.0003072	NA	11	2208.216	-4394.325	68.86300	1.111284e-15
30	-2.293086	-8.0e-06	0.0006266	-1e-07	-0.0006285	0.0004133	0.0002597	NA	0.0000657	-1e-07	0.0003026	-0.0001694	12	2206.324	-4388.520	74.66758	6.100666e-17
23	-2.292138	-8.1e-06	0.0005912	-1e-07	-0.0006038	0.0003798	-0.0003934	-1.4e-06	0.0001065	-1e-07	NA	-0.0001426	12	2202.054	-4379.981	83.20700	8.532288e-19
29	-2.293309	-8.1e-06	0.0006369	-1e-07	-0.0006437	0.0004255	0.0023796	-7.0e-07	NA	-1e-07	0.0003014	-0.0001575	12	2201.995	-4379.862	83.32553	8.041324e-19
15	-2.293107	-8.3e-06	0.0006330	-1e-07	-0.0006269	0.0004432	-0.0007647	-1.7e-06	0.0001292	-1e-07	0.0003042	NA	12	2201.708	-4379.289	83.89845	6.038371e-19
31	-2.293228	-8.2e-06	0.0006348	-1e-07	-0.0006332	0.0004243	-0.0003138	-1.2e-06	0.0001027	-1e-07	0.0003006	-0.0001152	13	2195.526	-4364.904	98.28412	4.540571e-22

kable(Subset.dredged.extinct.model, "html", caption = "BAMM Extinction Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM Extinction Male-bias Sexual Selection Dredge Table

	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	Sexual_selection_ppca	bioclim4:Sexual_selection_ppca	log(range.size.m2):Sexual_selection_ppca	NPP:Sexual_selection_ppca	PC1.LIG:Sexual_selection_ppca	residuals.PC1:Sexual_selection_ppca	df	logLik	AICc	delta	weight
0	-2.291100	-7.6e-06	0.0005471	-1e-07	-0.0005670	0.0004011	0.0009422	NA	NA	NA	NA	NA	8	2239.623	-4463.188	0.00000	9.982421e-01
16	-2.291695	-7.7e-06	0.0005655	-1e-07	-0.0005911	0.0003593	0.0008204	NA	NA	NA	NA	-0.0002658	9	2233.363	-4448.652	14.53554	6.964394e-04
8	-2.292345	-7.7e-06	0.0005965	-1e-07	-0.0006020	0.0004437	0.0009637	NA	NA	NA	0.0003152	NA	9	2233.113	-4448.152	15.03575	5.423296e-04
2	-2.290976	-7.6e-06	0.0005418	-1e-07	-0.0005524	0.0003885	-0.0024843	NA	0.0001231	NA	NA	NA	9	2233.054	-4448.035	15.15316	5.114101e-04
1	-2.291399	-7.9e-06	0.0005605	-1e-07	-0.0005809	0.0004077	0.0013441	-1.2e-06	NA	NA	NA	NA	9	2228.666	-4439.258	23.92997	6.351986e-06
24	-2.292817	-7.8e-06	0.0006107	-1e-07	-0.0006217	0.0004057	0.0008568	NA	NA	NA	0.0003034	-0.0002315	10	2226.838	-4433.586	29.60133	3.727258e-07
18	-2.291559	-7.6e-06	0.0005602	-1e-07	-0.0005786	0.0003533	-0.0016797	NA	0.0000902	NA	NA	-0.0002452	10	2226.797	-4433.504	29.68333	3.577527e-07
4	-2.291653	-7.9e-06	0.0005740	-1e-07	-0.0005907	0.0004065	0.0023362	NA	NA	-2e-07	NA	NA	9	2225.687	-4433.301	29.88723	3.230772e-07
10	-2.292220	-7.7e-06	0.0005912	-1e-07	-0.0005877	0.0004314	-0.0023393	NA	0.0001187	NA	0.0003140	NA	10	2226.542	-4432.994	30.19380	2.771635e-07
17	-2.291756	-7.8e-06	0.0005687	-1e-07	-0.0005941	0.0003635	0.0009507	-4.0e-07	NA	NA	NA	-0.0002518	10	2222.457	-4424.824	38.36346	4.663555e-09
3	-2.291393	-7.9e-06	0.0005605	-1e-07	-0.0005682	0.0003937	-0.0032200	-1.9e-06	0.0001720	NA	NA	NA	10	2222.169	-4424.249	38.93904	3.497280e-09
9	-2.292576	-7.9e-06	0.0006068	-1e-07	-0.0006129	0.0004487	0.0012922	-1.0e-06	NA	NA	0.0003117	NA	10	2222.151	-4424.212	38.97548	3.434136e-09
20	-2.292081	-7.9e-06	0.0005863	-1e-07	-0.0006079	0.0003718	0.0020931	NA	NA	-1e-07	NA	-0.0002167	10	2219.415	-4418.741	44.44700	2.226858e-10
26	-2.292680	-7.7e-06	0.0006054	-1e-07	-0.0006090	0.0003997	-0.0016523	NA	0.0000905	NA	0.0003035	-0.0002108	11	2220.273	-4418.437	44.75033	1.913486e-10
12	-2.292870	-8.0e-06	0.0006222	-1e-07	-0.0006246	0.0004485	0.0023252	NA	NA	-1e-07	0.0003115	NA	10	2219.171	-4418.253	44.93462	1.745048e-10
6	-2.291537	-7.8e-06	0.0005690	-1e-07	-0.0005789	0.0003970	-0.0002205	NA	0.0000895	-1e-07	NA	NA	10	2219.115	-4418.140	45.04770	1.649120e-10
19	-2.291677	-7.8e-06	0.0005670	-1e-07	-0.0005823	0.0003628	-0.0022333	-1.1e-06	0.0001231	NA	NA	-0.0001991	11	2215.992	-4409.877	53.31107	2.647794e-12
25	-2.292858	-7.8e-06	0.0006128	-1e-07	-0.0006237	0.0004086	0.0009470	-3.0e-07	NA	NA	0.0003030	-0.0002218	11	2215.932	-4409.756	53.43142	2.493164e-12
5	-2.292014	-8.2e-06	0.0005902	-1e-07	-0.0006074	0.0004143	0.0028381	-1.4e-06	NA	-2e-07	NA	NA	10	2214.737	-4409.384	53.80399	2.069413e-12
11	-2.292555	-8.0e-06	0.0006063	-1e-07	-0.0006006	0.0004350	-0.0029759	-1.7e-06	0.0001609	NA	0.0003077	NA	11	2215.650	-4409.191	53.99633	1.879667e-12
28	-2.293197	-8.0e-06	0.0006311	-1e-07	-0.0006382	0.0004180	0.0021203	NA	NA	-1e-07	0.0003025	-0.0001829	11	2212.889	-4403.671	59.51672	1.189444e-13
22	-2.291971	-7.9e-06	0.0005818	-1e-07	-0.0005983	0.0003671	0.0002469	NA	0.0000652	-1e-07	NA	-0.0002034	11	2212.849	-4403.591	59.59667	1.142835e-13
14	-2.292756	-7.9e-06	0.0006172	-1e-07	-0.0006133	0.0004393	-0.0001273	NA	0.0000859	-1e-07	0.0003107	NA	11	2212.598	-4403.089	60.09851	8.892208e-14
21	-2.292218	-8.1e-06	0.0005934	-1e-07	-0.0006146	0.0003809	0.0024002	-8.0e-07	NA	-1e-07	NA	-0.0001865	11	2208.522	-4394.936	68.25187	1.508452e-15
27	-2.292778	-7.9e-06	0.0006110	-1e-07	-0.0006121	0.0004077	-0.0021330	-9.0e-07	0.0001191	NA	0.0003020	-0.0001710	12	2209.466	-4394.805	68.38296	1.412746e-15
7	-2.291971	-8.2e-06	0.0005884	-1e-07	-0.0005954	0.0004025	-0.0009542	-2.0e-06	0.0001395	-1e-07	NA	NA	11	2208.232	-4394.357	68.83056	1.129455e-15
13	-2.293157	-8.2e-06	0.0006351	-1e-07	-0.0006382	0.0004545	0.0027477	-1.2e-06	NA	-2e-07	0.0003072	NA	11	2208.216	-4394.325	68.86300	1.111284e-15
30	-2.293086	-8.0e-06	0.0006266	-1e-07	-0.0006285	0.0004133	0.0002597	NA	0.0000657	-1e-07	0.0003026	-0.0001694	12	2206.324	-4388.520	74.66758	6.100666e-17
23	-2.292138	-8.1e-06	0.0005912	-1e-07	-0.0006038	0.0003798	-0.0003934	-1.4e-06	0.0001065	-1e-07	NA	-0.0001426	12	2202.054	-4379.981	83.20700	8.532288e-19
29	-2.293309	-8.1e-06	0.0006369	-1e-07	-0.0006437	0.0004255	0.0023796	-7.0e-07	NA	-1e-07	0.0003014	-0.0001575	12	2201.995	-4379.862	83.32553	8.041324e-19
15	-2.293107	-8.3e-06	0.0006330	-1e-07	-0.0006269	0.0004432	-0.0007647	-1.7e-06	0.0001292	-1e-07	0.0003042	NA	12	2201.708	-4379.289	83.89845	6.038371e-19
31	-2.293228	-8.2e-06	0.0006348	-1e-07	-0.0006332	0.0004243	-0.0003138	-1.2e-06	0.0001027	-1e-07	0.0003006	-0.0001152	13	2195.526	-4364.904	98.28412	4.540571e-22

#Run model for DR
Subset.MCC.DR.top <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.DR.top, 'data/Subset.MCC.DR.top.rds')

#Run model for ND
Subset.MCC.ND.top <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.ND.top, 'data/Subset.MCC.ND.top.rds')

#Run the 100 models for DR and ND using the best model:
Subset.data.noMCC <- SS.subset %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP)

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Subset.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC, 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.DR.model.data <- lapply(Subset.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Subset.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(SS.subset$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Subset.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(0.9711, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Subset.DR.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.DR.pgls.models, "data/Subset.DR.pgls.models.rds")

#Now for Node Density:
Subset.ND.model.data <- lapply(nd.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC, 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.ND.model.data <- lapply(Subset.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Subset.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Subset.ND.model.data, y = Subset.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.ND.pgls.models, "data/Subset.ND.pgls.models.rds")

#BAMM Top Models 

Subset.MCC.Lambda.top <- gls(mean.lambda ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.lambda),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Lambda.top, 'data/Subset.MCC.Lambda.top.rds')

Subset.MCC.Mu.top <- gls(mean.mu ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ sqrt(var.mu),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Mu.top, 'data/Subset.MCC.Mu.top.rds')

#Now for BAMM 100 models

Subset.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.BAMM.model.data <- lapply(Subset.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Subset.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(mean.lambda ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.lambda),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.lambda.pgls.models, "data/Subset.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Subset.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(mean.mu ~ Sexual_selection_ppca
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ sqrt(var.mu),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.mu.pgls.models, "data/Subset.BAMM.mu.pgls.models.rds")

The phylogenetic signal from the DR model is ~0.97. The other models had phylogenetic signals ~ 1, as such they were fixed at 1 to avoid problems of convergence.

Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')
Subset.MCC.DR.top[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()

• DR lambda: 0.9711

Subset.DR.pgls.models <- readRDS("data/Subset.DR.pgls.models.rds")
Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')

Subset.cols <- brewer.pal(n = 7, name = "Dark2")[-6]

Subset.DR.pgls.summary <- lapply(Subset.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.DR.summary <- data.frame(Subset.MCC.DR.top$coefficients, confint(Subset.MCC.DR.top)) %>% tibble::rownames_to_column()

Subset.DR.pgls.summary <- bind_rows(Subset.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `Sexual_selection_ppca` = "Sexual Selection"
                    )

Subset.DR.pgls.summary$Parameter = factor(Subset.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.DR.summary$Parameter = factor(Subset.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.DR.plot <-Subset.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.DR.plot <- Subset.DR.plot + geom_errorbarh(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS DR")

####################################
#For ND

Subset.ND.pgls.models <- readRDS("data/Subset.ND.pgls.models.rds")
Subset.MCC.ND.top <- readRDS('data/Subset.MCC.ND.top.rds')

Subset.ND.pgls.summary <- lapply(Subset.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.ND.summary <- data.frame(Subset.MCC.ND.top$coefficients, confint(Subset.MCC.ND.top)) %>% tibble::rownames_to_column()

Subset.ND.pgls.summary <- bind_rows(Subset.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.ND.pgls.summary$Parameter = factor(Subset.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.ND.summary$Parameter = factor(Subset.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.ND.plot <-Subset.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.ND.plot <- Subset.ND.plot + geom_errorbarh(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS ND")

######################

#For Lambda
Subset.BAMM.lambda.pgls.models <- readRDS("data/Subset.BAMM.lambda.pgls.models.rds")
Subset.MCC.Lambda.top <- readRDS('data/Subset.MCC.Lambda.top.rds')

Subset.Lambda.pgls.summary <- lapply(Subset.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.Lambda.summary <- data.frame(Subset.MCC.Lambda.top$coefficients, confint(Subset.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Subset.Lambda.pgls.summary <- bind_rows(Subset.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.Lambda.pgls.summary$Parameter = factor(Subset.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Lambda.summary$Parameter = factor(Subset.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.Lambda.pgls.summary.RO <- dcast(Subset.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter)

## Using UCI as value column: use value.var to override.

Subset.Lambda.pgls.summary.RO$Estimate <- NULL
Subset.Lambda.pgls.summary.RO <- sapply(Subset.Lambda.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
Subset.Lambda.pgls.summary.RO <-melt(Subset.Lambda.pgls.summary.RO) %>% na.omit()
Subset.Lambda.pgls.summary.RO$Var1 <- NULL
colnames(Subset.Lambda.pgls.summary.RO) <- c("Parameter", "Estimate")


Subset.Lambda.plot <-Subset.Lambda.pgls.summary.RO %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Lambda.plot <- Subset.Lambda.plot + geom_errorbarh(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Speciaton")

######################

#For Mu

Subset.BAMM.Mu.pgls.models <- readRDS("data/Subset.BAMM.mu.pgls.models.rds")
Subset.MCC.Mu.top <- readRDS('data/Subset.MCC.Mu.top.rds')

Subset.Mu.pgls.summary <- lapply(Subset.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.Mu.summary <- data.frame(Subset.MCC.Mu.top$coefficients, confint(Subset.MCC.Mu.top)) %>% tibble::rownames_to_column()

Subset.Mu.pgls.summary <- bind_rows(Subset.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.Mu.pgls.summary$Parameter = factor(Subset.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Mu.summary$Parameter = factor(Subset.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))


Subset.Mu.pgls.summary.RO <- dcast(Subset.Mu.pgls.summary %>% filter(Parameter != "(Intercept)"), Estimate ~ Parameter)

## Using UCI as value column: use value.var to override.

Subset.Mu.pgls.summary.RO$Estimate <- NULL
Subset.Mu.pgls.summary.RO <- sapply(Subset.Mu.pgls.summary.RO, function(x) {
  remove_outliers(x, na.rm = T)})
Subset.Mu.pgls.summary.RO <-melt(Subset.Mu.pgls.summary.RO) %>% na.omit()
Subset.Mu.pgls.summary.RO$Var1 <- NULL
colnames(Subset.Mu.pgls.summary.RO) <- c("Parameter", "Estimate")

Subset.Mu.plot <-Subset.Mu.pgls.summary.RO %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Mu.plot <- Subset.Mu.plot + geom_errorbarh(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Extinction")

#########################################

grid.arrange(symmetrise_scale(Subset.DR.plot, "x"),
             symmetrise_scale(Subset.ND.plot, "x"),
             symmetrise_scale(Subset.Lambda.plot, "x"), 
             symmetrise_scale(Subset.Mu.plot, "x"), 
             nrow = 1)

Figure S12: Measures of sexual selection a phylogenetic pca of sexual dimorphism, social polygyny and paternal care provide some evidence to suggest sexual selection predicts extinction rate (DR and ND models). This figure is seen in Figure 1 within the manuscript

Table S14: MCC model estimates from the above model are presented here as numberical values with 95 % CIs.

Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC DR Estimates from Male-bias SS")

MCC DR Estimates from Male-bias SS

Parameter	Estimate	LCI	UCI
Sexual_selection_ppca	0.03887	0.009472	0.06827
log(range.size.m2)	-0.01131	-0.02055	-0.002058
bioclim4	2.448e-05	-6.115e-05	0.0001101
residuals.PC1	0.0004409	-0.00859	0.009472
PC1.LIG	0.006881	-0.001666	0.01543
NPP	-7.335e-07	-3.45e-06	1.984e-06

Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC ND Estimates from Male-bias SS")

MCC ND Estimates from Male-bias SS

Parameter	Estimate	LCI	UCI
Sexual_selection_ppca	0.0004381	-0.0004178	0.001294
log(range.size.m2)	-0.00018	-0.0004838	0.0001238
bioclim4	-4.586e-07	-3.601e-06	2.683e-06
residuals.PC1	7.875e-05	-0.0002533	0.0004108
PC1.LIG	0.0001073	-0.0001577	0.0003723
NPP	-1.723e-08	-1.229e-07	8.84e-08

Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Speciation Estimates from Male-bias SS")

MCC BAMM Speciation Estimates from Male-bias SS

Parameter	Estimate	LCI	UCI
Sexual_selection_ppca	0.0009422	-0.002845	0.004729
log(range.size.m2)	0.0005471	-0.0006171	0.001711
bioclim4	-7.604e-06	-1.939e-05	4.177e-06
residuals.PC1	0.0004011	-0.0008461	0.001648
PC1.LIG	-0.000567	-0.001625	0.0004906
NPP	-5.98e-08	-4.751e-07	3.555e-07

Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)") %>% pander(caption = "MCC BAMM Extinction Estimates from Male-bias SS")

MCC BAMM Extinction Estimates from Male-bias SS

Parameter	Estimate	LCI	UCI
Sexual_selection_ppca	0.003443	-0.004091	0.01098
log(range.size.m2)	0.0004388	-0.002251	0.003129
bioclim4	-1.361e-05	-4.164e-05	1.442e-05
residuals.PC1	0.000148	-0.002828	0.003124
PC1.LIG	-0.0009247	-0.003254	0.001405
NPP	-2.112e-07	-1.169e-06	7.471e-07

Table S15: Highest posterior density intervals are calculated for the above figure (ie; models using Male-bias sexual selection measures). using this measure the models using \(\lambda_{DR}\) have model estimates that fall on the positive side 95 % of the time in 100 trees. For \(\lambda_{ND}\), there is a positive skew, however the 95 % HPD interval overlaps zero. \(\lambda_{BAMM}\) also shows a lesser positive skew. The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -). These intervals are calculated in the same way as in Table S7 and Table S9.

Subset.hpd.DR.top <- list()
Subset.DR.pgls.summary <- na.omit(Subset.DR.pgls.summary)
for(x in unique(Subset.DR.pgls.summary$Parameter)) {
Subset.hpd.DR.top[[x]] = hdi(Subset.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.DR.top <- bind_rows(Subset.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Subset.hpd.DR.top, 'data/Subset.hpd.DR.top.rds')

#For ND
Subset.hpd.ND.top <- list()
Subset.ND.pgls.summary <- na.omit(Subset.ND.pgls.summary)
for(x in unique(Subset.ND.pgls.summary$Parameter)) {
Subset.hpd.ND.top[[x]] = hdi(Subset.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.ND.top <- bind_rows(Subset.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.ND.top, 'data/Subset.hpd.ND.top.rds')

###############
Subset.hpd.Lambda.top <- list()
Subset.Lambda.pgls.summary <- na.omit(Subset.Lambda.pgls.summary)
for(x in unique(Subset.Lambda.pgls.summary$Parameter)) {
Subset.hpd.Lambda.top[[x]] = hdi(Subset.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Lambda.top <- bind_rows(Subset.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Subset.hpd.Lambda.top, 'data/Subset.hpd.Lambda.top.rds')

################
Subset.hpd.Mu.top <- list()
Subset.Mu.pgls.summary <- na.omit(Subset.Mu.pgls.summary)
for(x in unique(Subset.Mu.pgls.summary$Parameter)) {
Subset.hpd.Mu.top[[x]] = hdi(Subset.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Mu.top <- bind_rows(Subset.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.Mu.top, 'data/Subset.hpd.Mu.top.rds')

Subset.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Subset DR HPD Interval")

Subset DR HPD Interval

	Sexual_selection_ppca	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	0.00911	-0.0158	-9.31e-05	0.000313	-0.0036	-2.96e-06
Upper	0.0608	-0.000501	2.76e-05	0.0128	0.00894	1.19e-06

Subset.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Subset ND HPD Interval")

Subset ND HPD Interval

	Sexual_selection_ppca	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.00031	-0.000426	-3.63e-06	-3.05e-05	-0.000249	-9.23e-08
Upper	0.00154	0.000166	1.64e-06	0.000487	0.000249	6.61e-08

Subset.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Speciation HPD Interval")

Subset BAMM Speciation HPD Interval

	Sexual_selection_ppca	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.013	-0.0126	-9.69e-05	-0.00745	-0.00989	-2.09e-06
Upper	0.0309	0.0166	5.42e-05	0.0187	0.0047	1.99e-06

Subset.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Extinction HPD Interval")

Subset BAMM Extinction HPD Interval

	Sexual_selection_ppca	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0516	-0.0222	-0.000423	-0.0328	-0.0186	-3.14e-06
Upper	0.0576	0.0306	0.000171	0.0207	0.0158	2.78e-06

Phylogenetic Path Analysis

We undertook a phylogenetic path analysis using a model of causal pathways a priori. To do this we use the phylopath function. Our path analysis is undertaken on the dataset subset to species with measures of male bias sexual selection. Additionallty, the path analysis is restricted to using one measure of speciation (\(\lambda_{DR}\)) and the MCC tree. Our reasoning for the causal links are described breifly below:

1. \(\lambda_{DR}\) dependent upon Sexual Dichromatism: Sexual dichromatism would directly impact speciation if there is an association with the evolvability of chromatism and speciation, this may more rapidly lead to niche divergence or reproductive isolation.
2. \(\lambda_{DR}\) dependent upon Sexual Selection: This was a major hypothesis in our study (see introduction).
   
3. \(\lambda_{DR}\) dependent upon Temperature Seasonality: More variable environments are broadly expected to impact the resources, niches and fitness of species; thus there is an expectation that more variable environments may impact speciation rate.
   
4. \(\lambda_{DR}\) dependent upon Range Size: As previously seen in the PGLS models range size shows a negative correlation with speciation rate; expected reasons for this causal link include lower dispersal ability in smaller ranges correlating with lower gene flow and greater speciation. It may also reflect niche specialisation and thus in line with the previous point increase speciation and reduce gene flow.
   
5. Sexual Dichromatism dependent upon Sexual Selection: Higher levels of sexual selection will promote showy males and drab females (sexual dichromatism).
   
6. Sexual Dichromatism dependent upon Temperature Seasonality: Ecological aspects may facilitate sexual dimorphism as niche partitioning between the sexes causes sex-limited traits and colouration that are independent of mate choice or male-male competition.
   
7. Sexual Selection dependent upon Temperature Seasonality: Greater environmental variability/stress may impact the usefulness (fitness effects) of sexual selection (see introduction).
   
8. Range Size dependent upon Sexual Selection: Sexual selection may often depend on local adaptations and the alignment of sexually selected traits with locally valuable traits, range size may limit this leading to a negative correlation. Alternatively if sexual selection has added fitness and evolvability benefits there may be increased range sizes in sexually selected species.
   
9. Range Size dependent upon Temperature Seasonality: Anticipated to be more of a correlation than causal link. Yet a speecies range size may be dependent upon how variable their environment is and the availability to local resources.

#Set rownames to match tree
rownames(SS.subset) <- SS.subset$TipLabel
pruned.MCC.Subset.tree <- readRDS('data/pruned.MCC.Subset.tree.rds')
SS.subset$log.range.size <- log(SS.subset$range.size.m2)

SS.subset2 <- SS.subset %>% rename(
  DR = MCC.DR,
  SD = SDi,
  TS = bioclim4,
  RS = log.range.size,
  SS = Sexual_selection_ppca
)

models.Subset <- define_model_set(
  one = c(DR ~ SD, 
          DR ~ SS,
          DR ~ TS,
          DR ~ RS,
          SD ~ SS,
          SD ~ TS,
          SS ~ TS,
          RS ~ SS,
          RS ~ TS)
)

result <- phylo_path(models.Subset, data = SS.subset2, tree = pruned.MCC.Subset.tree, model = 'lambda')
  
best_model <- best(result)
set.seed(1)
path.plot <- plot(x = best_model,
     type = "color",
     algorithm = 'gem',
     manual_layout = NULL,
     curvature = 0.1,
     colors = c("#b2182b", "#2166ac"),
     show.legend = F)

#Inspect the coefficients with their SE

coef_plot(best_model, error_bar = "se", reverse_order = TRUE) + 
  ggplot2::coord_flip()+
  ggplot2::theme_minimal()

# pdf("Figures/Path_Plot.pdf", width=8, height=8)
# path.plot
# dev.off()

#path.plot

Figure S13: Path analysis standardized regression coefficients vary across relationships. Paths with error bars not everlapping zero are presented with an asterix in Figure 3 within the manuscript.

Additional Figures and Tables

BAMM.df <- readRDS('data/BAMM.df.rds')
BAMM.Rates.100Trees <- do.call(rbind, BAMM.df)
BAMM.Rates.100Trees$mean.lambda <- BAMM.Rates.100Trees$mean.lambda
BAMM.Lambda.100Trees <- ggiraphExtra::summarySE(BAMM.Rates.100Trees, measurevar = "mean.lambda", groupvars = "TipLabel")

DR.list <- plyr::llply(es.list, function(x) {
  x %>% as.data.frame %>% rownames_to_column()
  })
DR.list <- do.call(rbind, DR.list)
DR.Lambda.100Trees <- ggiraphExtra::summarySE(DR.list, measurevar = ".", groupvars = "rowname")
DR.Lambda.100Trees <- DR.Lambda.100Trees %>% rename(DR.Lambda.100Trees, DR = ., TipLabel = rowname)

BAMM.Lambda.100Trees %>% ggplot(aes(x = mean.lambda, y = DR.Lambda.100Trees$DR))+
  geom_point()+
  geom_errorbarh(aes(xmin = mean.lambda - 1.96*sd, xmax = mean.lambda + 1.96*sd), size = 0.05)+
  geom_errorbar(aes(ymin = DR.Lambda.100Trees$DR - 1.96*DR.Lambda.100Trees$sd, ymax = DR.Lambda.100Trees$DR + 1.96*DR.Lambda.100Trees$sd), size = 0.05)+
  theme_minimal()

# mean(DR.Lambda.100Trees$sd)
# mean(BAMM.Lambda.100Trees$sd)

Figure S14:

MCC.BAMM.model.data <- readRDS('data/MCC.BAMM.model.data.rds')
restricted.data %>% ggplot(aes(x = MCC.DR, y = MCC.BAMM.model.data$mean.lambda))+
  geom_point()+
  theme_minimal()

Figure S15:

What if we run the subset df qwith SDi

R Session information

sessionInfo() %>% pander()

R version 3.3.1 (2016-06-21)

**Platform:** x86_64-apple-darwin13.4.0 (64-bit)

locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8

attached base packages: parallel, grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: bindrcpp(v.0.2.2), ggtree(v.1.6.11), EBImage(v.4.16.0), phylopath(v.1.0.1.9000), ggraph(v.1.0.1.9999), stargazer(v.5.2.2), tibble(v.1.3.4), RColorBrewer(v.1.1-2), kableExtra(v.0.7.0), phangorn(v.2.4.0), HDInterval(v.0.2.0), MuMIn(v.1.42.1), coda(v.0.19-1), car(v.3.0-2), carData(v.3.0-1), BAMMtools(v.2.1.6), dplyr(v.0.7.6), ggExtra(v.0.8), reshape2(v.1.4.3), purrr(v.0.2.4), caper(v.1.0.1), mvtnorm(v.1.0-6), MASS(v.7.3-50), ggridges(v.0.5.0), brms(v.2.4.0), Rcpp(v.0.12.18), phytools(v.0.6-44), maps(v.3.2.0), gridExtra(v.2.3), geiger(v.2.0.6), lme4(v.1.1-18-1), Matrix(v.1.2-14), repmis(v.0.5), diversitree(v.0.9-10), ape(v.5.1), mgcv(v.1.8-24), nlme(v.3.1-128), tidyr(v.0.7.2), knitr(v.1.20), pander(v.0.6.2) and ggplot2(v.3.0.0.9000)

loaded via a namespace (and not attached): readxl(v.1.0.0), backports(v.1.1.2), fastmatch(v.1.1-0), plyr(v.1.8.4), igraph(v.1.1.2), lazyeval(v.0.2.1), splines(v.3.3.1), crosstalk(v.1.0.0), rstantools(v.1.5.1), inline(v.0.3.15), digest(v.0.6.16), htmltools(v.0.3.6), tiff(v.0.1-5), viridis(v.0.5.1), rsconnect(v.0.8.8), gdata(v.2.18.0), magrittr(v.1.5), openxlsx(v.4.1.0), readr(v.1.1.1), matrixStats(v.0.54.0), R.utils(v.2.7.0), xts(v.0.11-0), jpeg(v.0.1-8), colorspace(v.1.3-2), ggrepel(v.0.8.0), rvest(v.0.3.2), haven(v.1.1.2), jsonlite(v.1.5), bindr(v.0.1.1), survival(v.2.39-4), zoo(v.1.8-3), glue(v.1.3.0), gtable(v.0.2.0), R.cache(v.0.13.0), rstan(v.2.17.3), BiocGenerics(v.0.20.0), abind(v.1.4-5), scales(v.1.0.0), miniUI(v.0.1.1.1), plotrix(v.3.7-3), viridisLite(v.0.3.0), xtable(v.1.8-2), units(v.0.5-1), foreign(v.0.8-66), subplex(v.1.4-1), deSolve(v.1.20), stats4(v.3.3.1), StanHeaders(v.2.17.2), DT(v.0.4), animation(v.2.5), htmlwidgets(v.1.2), httr(v.1.3.1), threejs(v.0.3.1), gplots(v.3.0.1), pkgconfig(v.2.0.2), loo(v.2.0.0), R.methodsS3(v.1.7.1), locfit(v.1.5-9.1), tidyselect(v.0.2.4), rlang(v.0.2.2), later(v.0.7.3), munsell(v.0.5.0), cellranger(v.1.1.0), tools(v.3.3.1), fftwtools(v.0.9-8), evaluate(v.0.10.1), stringr(v.1.2.0), yaml(v.2.2.0), tidygraph(v.1.1.0), zip(v.1.0.0), caTools(v.1.17.1), mime(v.0.5), R.oo(v.1.22.0), xml2(v.1.2.0), bayesplot(v.1.6.0), shinythemes(v.1.1.1), png(v.0.1-7), curl(v.3.2), clusterGeneration(v.1.3.4), tweenr(v.0.1.5), stringi(v.1.1.6), Brobdingnag(v.1.2-6), forcats(v.0.3.0), lattice(v.0.20-33), nloptr(v.1.0.4), markdown(v.0.8), shinyjs(v.1.0), msm(v.1.6.6), combinat(v.0.0-8), bridgesampling(v.0.5-2), data.table(v.1.11.4), bitops(v.1.0-6), httpuv(v.1.4.5), R6(v.2.2.2), promises(v.1.0.1), KernSmooth(v.2.23-15), rio(v.0.5.10), codetools(v.0.2-14), colourpicker(v.1.0), gtools(v.3.8.1), assertthat(v.0.2.0), rprojroot(v.1.3-2), withr(v.2.1.2), shinystan(v.2.4.0), mnormt(v.1.5-5), expm(v.0.999-2), hms(v.0.4.2), udunits2(v.0.13), quadprog(v.1.5-5), minqa(v.1.2.4), rmarkdown(v.1.10), ggforce(v.0.1.1), numDeriv(v.2016.8-1), scatterplot3d(v.0.3-41), shiny(v.1.1.0), base64enc(v.0.1-3) and dygraphs(v.1.1.1.6)

References

Armenta, J. K., P. O. Dunn, and L. A. Whittingham. 2008. Quantifying avian sexual dichromatism: A comparison of methods. Journal of Experimental Biology 211:2423. Journal Article.

BirdLife International and Handbook of the Birds of the World. 2017. Bird species distribution maps of the world. http://datazone.birdlife.org/species/requestdis.

Dale, J., C. J. Dey, K. Delhey, B. Kempenaers, and M. Valcu. 2015. The effects of life history and sexual selection on male and female plumage colouration. Nature 527:367–370. Journal Article.

Del Hoyo, J., A. Elliott, and D. Christie. 2011. Handbook of the birds of the world (Vols. 8-16). Book, Lynx Edicions 2003-2011.

Harvey, M. G., G. F. Seeholzer, B. T. Smith, D. L. Rabosky, A. M. Cuervo, and R. T. Brumfield. 2017. Positive association between population genetic differentiation and speciation rates in new world birds. Proceedings of the National Academy of Sciences 114:6328–6333.

Harvey Michael, G., L. Rabosky Daniel, and N. Cooper. 2017. Continuous traits and speciation rates: Alternatives to state‐dependent diversification models. Methods in Ecology and Evolution 9:984–993. Journal Article.

Jetz, W., G. H. Thomas, J. B. Joy, K. Hartmann, and A. O. Mooers. 2012. The global diversity of birds in space and time. Nature 491:444–448. Journal Article.

Quintero, I., and W. Jetz. 2018. Global elevational diversity and diversification of birds. Nature 555:246.

Schliep, K. P. 2010. Phangorn: Phylogenetic analysis in r. Bioinformatics 27:592–593.